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SUMMARY 


A new method has been developed for two and three-dimensional 
computations of viscous supersonic flows with embedded subsonic regions 
adjacent to solid boundaries. The approach employs a reduced form of the 
Navier-Stokes equations which allows solution as an initial-boundary value 
problem in space, using an efficient noniterative forward marching 
algorithm. Numerical instability associated with forward marching algorithms 
for f lows with embedded subsonic regions is avoided by approximation of the 
reduced form of the Navier-Stokes equations in the subsonic regions of the 
boundary layers. Supersonic and subsonic portions of the flow field are 
simultaneous calculated by a consistently split linearized block implicit 
computational algorithm. The results of computations for a series of test 
cases relevant to internal supersonic flow is presented and compared with 
data* Comparison between data and computation are in general excellent thus 
indicating that the computational technique has great promise as a tool for 
calculating supersonic flow with embedded subsonic regions. Finally, a 
User f 8 Manual is presented for the computer code used to perform the 
calculations . 


INTRODUCTION 


The interaction of a supersonic jet with an external stream has received 
considerable attention due to its critical role in determining both the 
thermodynamic and acoustic signatures from a vehicle producing a supersonic 
jet. The physics of supersonic jet flow is characterized by the formation of 
shock and expansion waves, the development of viscous mixing layers and the 
interaction of these phenomena (Figs. 1 and 2). As a consequence, a 
complicated multiple shock cell flow structure comes into existence in the 
jet. The detailed flow structure of a given jet including a shape of jet 
interface is determined primarily by the conditions of an ambient flow. At 
the jet interface, these two different flows form a turbulent mixing layer 
which strongly interacts with both shock and expansion waves. In the 
presence of interactions with waves, the turbulence structure in the mixing 
layer may be strongly affected. Both compressibility effects and pressure 
gradients in the jet can also influence the turbulence structure in the 
mixing layer. Under certain operational conditions, Mach disc(s) may occur 
and complicate the turbulent jet flow structure (Fig. 2). The mixing process 
that is present downstream of the disc has wake-like characteristics, and the 
entrainment process induces large streamline displacement which can 
appreciably alter pressure levels from those occurring in the inviscid 
limit. Because of the strong viscous/ inviscid interactions as discussed 
above, inviscid theoretical methods have been found to be inadequate for 
predicting jet/external flow interactions. Therefore, very costly 
experimental testing has been the primary means of determining the complex 
flow patterns associated with these flows. Hov.ever, in recent years, 
advances in computational fluid dynamics have demonstrated that the 
dependence on the experimental analysis of such a complex flow can be reduced 
by utilizing computational analysis. Unlike the experimental testing, 
computational analysis is not necessarily restricted by Reynolds number and 
other flow or testing conditions. Furthermore, experimental costs are 
steadily increasing, whereas the computational costs have been and are 
expected to decrease. Thus, the goal of the current work was to develop an 
efficient and reliable numerical approach for the prediction of the detailed 
flow structure of a supersonic jet exhausted into the external flow. 


2 



Early computational viscous solutions of the nozzle flow field region 
consisted of viscous-inviscid interaction techniques called patching methods 
(for example, Refs, 1-3). The techniques consist of dividing the various 
viscous and inviscid regions of the flow according to the physical nature and 
mathematical behavior of the equations applicable to each region. Each of 
these regions was then analyzed independently and was iteratively coupled 
using appropriate matching conditions consistent with weak^interac tion 
theory. Although these solutions gave reasonable results for specific data 
sets, the required matching limited the usefulness of this method as a 
predictive technique. It is quite natural to consider solving the full 
time-dependent, compressible Navier-Stokes equations as an alternative 
approach for such a complex flow (for example, Refs. 4-5). However, the 
widely disparate length scales and flow characteristics in the various 
regions involved could lead to perhaps unacceptable computer time and storage 
requirements in achieving the results of adequate resolution, especially in 
three dimensions. In this study, it was decided that such a procedure should 
be used only if no suitable alternative exists, A compromise approach would 
possess the general three-dimensional, viscous nature of the Navier-Stokes 
equations, but would take advantage of realistic physical approximations to 
limit the computer running time and storage requirements associated with the 
solution of the complete three-dimensional Navier-Stokes equations. 

During the past two decades much effort has been expended in developing 
such numerical procedures which can be used as an alternative to solving the 
full Navier-Stokes equations for certain classes of problems. These 
procedures treat a reduced form of the steady state Navier-Stokes equations, 
often referred to as the parabolized Navier-Stokes equations, as an initial 
boundary value problem that can be solved by spatial forward marching. The 
ability to forward march the governing equations from an initial streamwise 
location to some desired downstream location rather than perform a 
simultaneous solution of the governing equations at all streamwise locations, 
as is required for the solution of the full Navier-Stokes equations, results 
in a considerable reduction of computational time. Although the amount of 
reduction will depend on the problem considered, the efficiency of the 
solution procedures and numerous other variables, this reduction has been the 
primary motivation for the development of marching procedures. 
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To devise a set of governing equations suitable for the spatial forward 
marching of supersonic flows, three steps must be taken. First, a nominally 
primary flow direction must be identified. Second, a coordinate system must 
be constructed with one of its coordinate directions aligned with the primary 
flow direction. Third, all diffusion in the primary flow direction must be 
neglected. These steps when applied to the steady state Navier-Stokes 
equations produce a set of governing equations which is well posed for the 
spatial forward marching of supersonic flows (e.g. Ref. 6). The introduction 
of no-slip surfaces into a supersonic flow usually results in the formation 
of embedded subsonic regions adjacent to these surfaces. Also, under certain 
operating conditions of the vehicle supersonic flow coexists with the 
subsonic flow as can be found in the case of supersonic jet exhaust into the 
subsonic ambient flow. When the set of reduced equations is forward marched 
with embedded or coflowing subsonic regions, the governing equations are not 
well posed and hence the solution procedure may become unstable. This 
instability, which is often referred to as the branching phenomenon, has been 
the subject of much research (e.g. Refs. 7-8) and the techniques used to 
suppress this instability are a convenient way to differentiate between 
procedures for solving the reduced form of the Navier-S tokes equations for 
supersonic flow with subsonic regions. In one of the earliest works in this 
area Garvine (Ref. 7) demonstrated (for a model problem) the existence of 
exponentially growing (divergent) terms in the solution of the reduced form 
of the Navier-Stokes equations when applied to a viscous (inviscid supersonic 
flow interacting with a viscous boundary layer) interaction problem. The 
author concluded that for this problem the reduced form of the Navier-Stokes 
equtions was improperly set as an initial value problem, because the 
interaction dynamics contained upstream ’'elliptic' 1 influence. In the model 
problem if the upstream conditions are not precisely set to cause the 
divergent terras to be multiplied by zero, the exponentially growing terms 
will cause the streamwise pressure gradient terms to grow exponentially large 
resulting in unrestrained acceleration or deceleration of the flow. In 
general it is not possible to pick the upstream conditions to negate the 
exponentially growing modes, hence several investigators have attempted to 
suppress the unstable (or branching) behavior by further modification of the 
reduced form of the Navier-Stokes equations for supersonic flows with 
embedded subsonic regions. 
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Much of the early work on the solution of the reduced form of the 
Navier-Stokes equations is based on the work of Rudman and Rubin (Ref. 8). 
Rudman and Rubin solved the equations for the hypersonic flow over slender 
bodies with sharp leading edges. Based on an order of magnitude analysis 
they demonstrated that for this class of problems the streamwise pressure 
gradient term was negligible when compared with the inertia and viscous terms 
of the streamwise momentum equation* Neglecting the streamwise pressure 
gradient term in this equation prevented the branching behavior in their 
calculation. Although this approach does yield a set of equations that is 
well posed for spatial forward marching, its assumption of negligible 
streamwise pressure gradient limits the cases which can be considered* In a 
later work Lubard and Helliwell (Ref. 9) proposed a method for preventing 
branching that involved explicit spatially lagged evaluation of the 
streamwise pressure gradient term. The above authors found that in addition 
to the frequently encountered problem of instability associated with 
exceeding some marching direction step size, a further instability was 
encountered when the step size was reduced below some limit. By examining 
the eigenvalues of a model set of equations (Ref. 12) they were able to 
develop a criterion for this minimum step size. Numerical experimentation 
with their computer code demonstrated reasonable co l l elation with their 
criterion. Several test cases using this method have been successfully run 
(mainly for cone flow cases) by the authors of Ref. 9 and others 
(Refs. 13-14) and in these cases evidently the restriction on the minimum 
marching step size was not a problem in allowing sufficiently accurate 
results to be obtained. However, the restriction on minimum marching step 
size is in principle not a desirable feature. Since it does prevent 
arbitrary mesh refinement, and thereby the assurance that an accurate unique 
solution has been obtained. In at least one case (Ref. 9) this minimum step 
size restriction prevented the successful running of a case. In a later 
technique developed by Vigneron, Rakich and Tannehill (Ref. 10) a variant of 
the technique of Lubard and Helliwell was used to prevent branching. In this 
particular variant the streamwise pressure gradient term was evaluated by an 
implicit backward difference in the supersonic portion of the flow. However, 
in the subsonic region only that portion of the streamwise pressure gradient 
term that could be included without causing branching was evaluated 
implicitly. The stability analysis through the investigation of the 
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eigenvalues of a model set of equations also produces a restriction on the 
minimum allowable step size. When that portion of the subsonic pressure 
gradient excluded from the implicit evaluation is evaluated explicitly by a 
lagged technique similar to Lubard and Helliwell, the investigators noted 
that the scheme became unstable. Thus, in order to achieve stability this 
technique neglected the explicit portion of the streamwise pressure gradient 
terra in the subsonic region. Schiff and Steger (Ref. 11) treat the subsonic 
streamwise pressure gradient terra by either a first- or second-order 
streamwise extrapolation technique in the subsonic regions. The first order 
technique is equivalent to setting the steamwise pressure gradient term equal 
to zero in the subsonic region while the second order technique is equivalent 
to the explicit evaluation of the streamwise pressure term (as was done by 
Lubard and Helliwell). As with the two previously discussed techniques, 
these authors also report a restriction on the minimum marching step size 
that they may take and still retain a stable calculation. On the other hand, 
Rubin and Lin (Ref. 15) have developed a global relaxation procedure for 
solving the reduced form of the Navier-Stokes equations. This technique was 
primarily developed for application to cases where upstream influence is 
strong. To obtain the upstream influence with the reduced form of the 
Navier-Stokes equations requires a global iteration or relaxation procedure. 
The above authors do this by approximating the streamwise pressure gradient 
term by a forward difference. When marching. the solution from the i th to 
the i+ ist station the pressure gradient terra is evaluated in terras of the 
pressure at the i+l st (current) and i+2 nc * (unknown) station (i+l st 
station is the implicit station; all other streamwise derivatives are 
backward differenced between i+1 and i) . Initially the pressure at the 
unknown i+2 n< * station is guessed; during subsequent global iterations the 
previously calculated value is used. Global iteration of the governing 
equations is continued until the solution converges. Rubin and Lin report 
that convergence is typically obtained in five to ten iterations for cases 
with small streamwise pressure gradients (cases run to date have been limited 
to flow over cones). The authors also report that there is no minimum 
marching step size requirement with their approach. 

The supersonic jet exhausted from a nozzle is essentially turbulent and 
has associated with it large streamwise pressure gradients, e.g., turbulent 
mixing layer interacting with multiple shock wave-cells. It is expected that 
one would desire to take a small streamwise marching step to accurately 
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resolve the phenomenon, for example, the near-field plume-mixing 
interactions . 

As reviewed above, the methods of Refs. 8-11 all make an attempt to 
consider the effects of streamwise pressure gradient in the embedded subsonic 
regions. However, they give only an approximate treatment to this possibly 
dominant term, and all of those methods have a minimum marching step size 
limitation which, in the cases of interest in this study, may not either 
allow for an accurate or in some cases even a minimally acceptable solution. 
Although there is no minimum marching step size requirement with the global 
relaxation procedure developed by Rubin and Lin (Ref. 15), the computation 
time would increase considerably in more complex flow situations. If the 
number of iterations were to increse much beyond the five to ten iterations 
it took to converge a small streamwise pressure gradient case, the 
computation time would then become comparable to that required for a solution 
of the full Navier-Stokes equations. 

Hence, a noniterative approach is sought with a consequent reduction in 
computer cost relative to either the global iteration approach to solving the 
reduced Navier-Stokes equations or solution of the full Navier-Stokes 
equation. Further as a prerequisite it is required that there exist no 
numerical limitation on the minimum marching step and it is desired to keep 
to a minimum any approximation to the streamwise pressure gradient term. 

Such an approach which meets these requirements has been developed by 
the present authors as reported elsewhere (Ref. 16). A detailed explanation 
of the techniques will be presented in the next analysis section. Unlike 
other noniterative methods, the streamwise pressure gradient terra is 
maintained and evaluated by an implicit backward difference in the supersonic 
as well as the subsonic region. However, approximations were made on both 
the normal (to the wall usually) momentum and continuity equation in the 
subsonic region(s) to suppress the branching behavior of solution. In the 
(thin) subsonic flow region, the normal momentum equation is approximated by 
the conventional boundary layer form of that equation, and the continuity 
equation is approximated by integrating the continuity equation from the wall 
to an arbitrary point in the subsonic portion of the boundary layer. This 
approximation of continuity equation is a major difference between the 
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present technique and other noniterative methods. As discussed in Ref. 16, 
numerical experimentation with various cases indicates no restriction on the 
minimum step size. Thus, an accurate resolution of the flow fields has been 
successfully obtained without being limited by the requirement of minimum 
step size. Since no approximations were made on the streamwise pressure 
gradient term, the predictions for the flow field under the condition of 
non-negligible streamwise pressure gradient agreed very well with the data as 
discussed in Ref. 16. 

Thus, application of this approach to the analysis of steady 
three-dimensional supersonic jet interacting with an ambient flow was a major 
concern in this study. The flow field calculation by the existing computer 
program, called PEPSIS, was limited to a computational domain without any 
embedded solid objects except for the wall boundary. Thus, the first task of 
this study was to extend a capability of PEPSIS so that the flow field could 
be analyzed in a computational domain with one or more embedded solid bodies 
such as nozzle walls. The extended PEPSIS computer program, which is called 
PEPSIN, is the same as PEPSIS as far as its basic assumptions and 
approximations are concerned. One of the basic assumptions for PEPSIS 
approach was the thinness of subsonic layer within a boundary layer. As a 
result of this assumption, approximations were made on the governing 
equations within this thin, subsonic layer based on the fundamental 
properties of a boundary layer. However, analysis of supersonic jet 
interactions under certain low conditions is expected to encounter subsonic 
regions of a different nature. One of such subsonic regions may be present 
when a Mach disc appears in the jet. Another type of subsonic region may be 
encountered when a supersonic jet exhausts into an ambient subsonic stream. 
First of all, these two kinds of subsonic regions may not be necessarily 
thin. Furthermore, since these subsonic regions are not located in the 
vicinity of a wall boundary, approximations based on boundary layer 
properties become questionable. Therefore, it was decided to exclude such 
subsonic regions from consideration in the present study. Those subsonic 
regions require further separate investigations before they can be 
incorporated in the PEPSIN approach. Thus, as will be discussed later, test 
cases considered in the present case were limited to the analysis of a 
supersonic jet interacting with supersonic external stream. 
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Two other investigators have utilized forward marching techniques to 
solve the jet/external flow interactions. First, Dash and co-workers 
(Refs. 17 and 18) developed an explicit spatial marching method to predict 
both the inviscid plume structure and the viscous mixing layer. Their 
approach employs a multiple-domain solution algorithm to solve a reduced form 
of the Navier-Stokes equations. Their approach has been primarily applied to 
the two-dimensional axisymmetric cases, although they presented some results 
for nonaxisymmetric three-dimensional cases (rectangular cross sections of 
various aspect ratios) in a recent paper (Ref. 19). However, the marching 
step size in their procedure is limited by numerical stability conditions of 
explicit schemes. On the other hand, Vatsa and his co-workers (Ref. 18) 
developed a method to solve the problem of a slightly underexpanded 
supersonic jet in a subsonic co-flowing mainstream. First, they solved the 
inviscid problem which was formulated in terms of linearized perturbation 
potentials to estimate the streamline curvature. Then they obtained 
governing equations of parabolic type through the use of an approxiarate 
inviscid flow intrinsic coordinate system. Their solution based on a 
marching technique produces the class of weak interactions wherein the 
viscous interaction has a formally second-order effect. In this approach, 
prediction of streamline curvature, in general, may become as difficult as 
the viscous problem if the difference between jet and external flow 
conditions become large or if the problem is three-dimensional. A difficulty 
may arise also if the viscous effects become significant and have a strong 
influence on the secondary flow. 

The present study is concerned with the application of an implicit 
forward marching procedure for the reduced form of the Navier-Stokes 
equations (Ref. 16) to the three-dimensional analysis of supersonic jet flow 
interacting with a co-flowing stream. The approach is noniterative and is 
not limited by a minimum step size criterion. Furthermore, being implicit 
the maximum allowable step size is expected to exceed that of explicit 
methods by a significant amount. The remainder of this report will consist 
of (1) a discussion of the analysis used in the study, (2) a discussion of 
the solution of the governing equations, (3) the results of a series of test 
cases run to demonstrate the applicability of the analysis and to exercise 
and to validate the resulting computer code and (A) a User’s Manual for the 
computer code, termed PEPSIN. 
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LIST OF SYMBOLS 


A 

C H 

C P 

D 

D:D 

H 

L 

M 

P 

Pr 

Re 

S 

T 

U 

V 

W 

X 

$ 

h 

> 

n 


Square Matrix of coefficients 
Nondimensional heat transfer coefficient 
Specific heat 

Column vector whose elements are spatial differential 
operators 

Second invariant of the mean flow rate of deformation tensor 

Column vector associated with marching direction terms 

Linear differential operator 

Mach number 

Static pressure 

Prandtl number 

Reynolds number 

Source term 

Static temperature 

Streamwise velocity component 

Transverse velocity component 

Spanwise velocity component 

Distance from leading edge 

Velocity 

Damping coefficient 
Metric coefficient, enthalpy 
Mixing length 

Unit vector normal 
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LIST OF SYMBOLS 


Heat transfer 
Radius 

Velocity component 
Coordinate direction 
Distance from a surface 
Nondiraensional distance from a surface 


Greek Symbols 

Ratio of specific heats 
Boundary layer thickness 
von Karman constant 
Viscosity 
Density 

Viscous stress tensor 
Nabla operator 

Marching direction step size 


Subscripts 


^thd direction 

j th direction, jet flow condition 
Laminar 

Normal direction 
Wall 

External stream condition 
Jet flow condition 
Red uc ed 
Sonic line 

Turbulent, tangential direction 
Stagnation condition 
Streamwise direction 
Cross plane direction 
Cross plane direction 
Free stream condition 


Superscripts 

i ttl streamwise station 


Transpose 


II. ANALYSIS 


Governing Equations 

The steady state fluid dynamic conservation laws of mass, momentum and 
energy respectively can be written in nondimensional vector form as 

V/>\7=0 (1) 


v* (ov7)+ vp- Yn= o 

r Re 


( 2 ) 


and 


VW>— *•[£(•£♦ -£l)vt]-v- 


(r* V) 
Re 


■= O 


( 3 ) 


This form of Che governing equations, often referred to as the full 
Navier-Stokes equations, requires several auxiliary relations and models 
before these equations can be solved. In this study the stagnation enthalpy, 

ho, Is related to the static temperature, T, and the velocity V, through the 
relationship (assuming a calorlcally perfect gas) 

"v • V 

h o s C p T+— g— (4) 


while the temperature, T, pressure P, and density, p, are related by means of 
the calorlcally perfect gas equation of state 

P =ZZic p ^T 

r p 

The stress tensor, t, is represented by the relationship 

r =/i.(W + VV t )~Y/aV*V 

where the superscript T refers to the transpose of the tensor. The 

components of the velocity vector, are Interpreted as the mass weighted 
mean velocity components and p, P and T are the ensemble-averaged 


( 5 ) 


( 6 ) 
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density, pressure and temperature (Ref. 21). Hence, these equations can be 
applied to both laminar and turbulent flows if the effective viscosity, p, is 
interpreted as the sum of the laminar, P£, and turbulent, 
viscosities , i .e . , 


H- * H-i + / i T (7) 

It is assumed that the laminar viscosity can be computed from Sutherland's 
law and that the laminar and turbulent Prandtl numbers, Pr£ and Pr-p are 
constant. For this study an algebraic mixing length turbulence model of the 
form 


/jl t = Re pl m VdTd 


( 8 ) 


was used where Am is the algebraic mixing length and D:D is the second 
invariant of the mean flow rate of deformation tensor (Ref. 22). In this 
study the mixing length of McDonald and Camarata (Ref. 23) was used. 

l m = 0.09 S b tanh [<y /(0.09 S b )]^> (9) 


where 6^ is the local boundary layer thickness, < is the von Karman 
constant, y is the distance to the nearest wall, and is t * ie sublayer 
damping terra of van Driest (Ref. 24). 

To obtain what is often referred to as the reduced form or the 
'parabolized' form of the Nav ier-Stokes equations involves approximation of 
the diffusion terms (both stress and Fourier heat conduction) of Eqs . (2), 
(3) and (6). This approximation neglects all derivatives of the stress 
tensor and the Fourier heat conduction terras in a selected 'marching' or 
'streamwise* direction. In addition all streamwise derivatives of the 
velocity components of the stress tensor are neglected. 

Hence, the general reduced form of Eqs. (2) and (3) can be recast as 


V-(pW) + VP- 


(V'-rh 


-= o 




Re 


( 10 ) 



and 




— 1 ) VT 1 } 

- fv- (rV) l 

«• Re ' 

. Pr, 

Pr T J J J R 

*• Re -I 


L* 


( 11 ) 


where the subscript R refers to the approximated or reduced form of the noted 
term. 

The reduced form of the Navier-Stokes equations, Eqs . (1), (10) and (11) 
is the starting point for the discussion of the governing equations to be 
used for this study. As discussed in detail in Ref. (16), this set of 
equations is not well posed for solution by spatial forward marching when 
applied to the class of problems considered in this study, i.e., supersonic 
flow with embedded subsonic layers. In this investigation, the same 
strategy as in Ref. (16) was adopted. The flow is divided into supersonic 
and subsonic flow regions and different sets of governing equations are 
utilized in each region. The reduced form of the Navier-Stokes equations 
in the supersonic region(s) of the flow and what can be considered to be a 
model set of equations in the subsonic region(s) are utilized. The model set 
of equations used in the subsonic region(s) is obtained by starting with the 
reduced form of the Navier-Stokes equation and making appropriate physical 
approximations in this region such that the coupled system of equations in 
both supersonic and subsonic flow are stable when solved as an initial value 
problem in space. There are several important features of the subsonic model 
set of governing equations. First, no approximation is made to either the 
streamwise or the tangential momenta equations, and second, the terms 
tangential and normal refer to some adjacent solid surface. Hence, the full 
effect of the pressure gradient terras will be felt in both the streamwise and 
tangential direction of the subsonic portion of the flow. The assumption 
needed to modify the normal momentum and continuity equations in the subsonic 
regions is the relatively unrestr ic tive condition that the subsonic layer is 
thin relative to the characteristic transverse dimension of the flow device. 
For the case of an impermeable wall this leads to the condition that within 
the viscous subsonic layer the normal velocity component is negligible. The 
importance of the specification of the normal velocity is that a mechanism 
can now be established to prevent the growing 
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mode caused by the interaction between the subsonic and supersonic layers, 
i.e., the branching phenomenon. This approximation can be obtained by 
integrating the continuity equation from the wall to an arbitrary point in 
the subsonic portion of the boundary layer X s , Thus, 


»i h T/>"n| s * ~^T^< h 2 h 3/> w i ) + a^ (h i h nf w T>] dx * +h i h T /> w l w 


( 12 ) 


where the subscripts n and T refer to the crossflow direction normal and 
tangential to the walls, s refers to the evaluation at the arbitrary point in 
the subsonic region, and the subscript W refers to the evaluation at the 
wall. Restricting our attention to flows where the subsonic region is 
sufficiently thin allows the integral in Eq • (12) to be neglected and hence 
this equation can be approximated by 

h , h T/> w n| s = h i h r/ )W n| w (13) 

As mentioned above, for the case of an impermeable wall Eq. (13) further 
reduces to 


w 


n 


* O 


(14) 


As a consequence of the thin subsonic layer assumption, the normal (to the 
wall) momentum equation then can be expressed as a balance only between the 
normal pressure gradient and the centrifugal (curvature) forces in the 
subsonic layer. For example, in general orthogonal coordinates, 

Xi , X 2 , and X 3 with corresponding metric coefficients h\, h 2 and h 3 and 
velocity components wi, W 2 and W 3 this equation is expressed as 


dP 

dx r 


’ T w T / ah 3 dh n v n w, / a^ _ah J1 \ 

n = ( -o hr ( w 3 ^ w n dx J + ( h, ' w « ax n " w " ax, ' 


(15) 


where and X^ respectively refers to the appropriate cross-sectional 
direction normal and tangential to the wall (n and T have values of 2 or 3; 
direction 1 is the nominally streamwise direction). In summary, the new set 
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of governing equations consisting of the reduced form of the Navier-Stokes 
equations in the supersonic portion of the flow and the model set of 
equations in the subsonic regions of the flow (Eq. 13 instead of continuity 
equation and Eq. 13 for the normal momentum equation) has, on the bases of 
numerical experimentation (Ref, 16), been found to be well posed for solution 
by spatial forward marching for a wide range of practical problems. 

Initial and Boundary Conations 

To uniquely define the problem of interest it is necessary to specify 
both initial and boundary conditions. For a spatial forward marching 
procedure, the initial conditions refer to the set of conditions that must be 
specified at the initial marching station. Boundary conditions must be set 
on the boundaries of the cross-sectional marching plane. For the calculation 
of the supersonic jet flow interacting with the coflowing stream, information 
must exist at an initial plane such that a reasonable approximation to a 
complete set of initial data can be constructed. In its most pure form this 
would be an initial plane where experimental or computational data was 
available such that all the initial conditions were known and consistent with 
the constitutive relations. In general, automatic generation of initial 
profiles is not straightforward for the present type of flow. However, if 
the flow is two-dimensional, a limited automation is possible. Usually a 
limited amount of information is available, where, for instance, freestream 
conditions with a boundary layer thickness and a skin friction coefficient on 
the nozzle surface might be known. In this case, a theoretical boundary 
layer profile of the pertinent variables (velocity components, temperature, 
pressure, etc.) can often be derived and matched with the flow outside the 
boundary layer. Analysis of the characteristics of the supersonic 
three-dimensional Euler equation shows that there are five characteristics 
entering the upstream boundary of the computational domain. It is argued 
that the reduced form of the Navier-Stokes equations used here would not 
require more boundary conditions than the Euler equations and this certainly 
appears reasonable at very high Reynolds numbers far from the solid walls. 
Thus, it is argued that five conditions must be set on the inflow boundary. 

In this study those conditions are chosen as the three velocity components, 
the pressure and the temperature. It is to be emphasized that the initial 
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conditions must in some sense be consistent with the governing equations in a 
hyperbolic system. In supersonic flow computations inconsistencies, 
perturbations, etc. can persist far downstream. 

As used in this investigation the boundary conditions utilized on the 
boundaries of the cross-sectional plane can be divided into three categories: 
(1) wall conditions, (2) symmetry conditions and (3) external flow 
conditions. Analysis of the characteristics of the boundary layer equations 
shows that four conditions must be specified on walls. Arguing that in the 
limit as the present system approaches the wall it should approximate the 
boundary layer equations, the same four conditions are used here. For this 
study the no-slip conditions are used for the strearawise and tangential cross 
plane velocity components, i.e., 


W, = 0 


( 16 ) 


and 


W T = 0 ( 17 ) 

where again the subscript 1 refers to the streamwise direction and the 
subscript T refers to the cross plane tangential velocity direction. For the 
cross plane normal velocity component either the normal velocity or the 
normal mass flux are specified, i.e., 


w 


n 


w w 


( 18 ) 


or 


>n = P w w w 


( 19 ) 


where the subscript W refers to the wall value (specified), 
condition (the thermal condition) used is to either specify 


The fourth 
an adiabatic wall 
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or to specify the wall temperature (a cold or hot wall/* The condition can 
be specified respectively as 


n w • VT = 0 


( 20 ) 


or 


T 


= T 


w 


( 21 ) 


where in this case n w represents the unit vector normal to the wall. In 
addition, a fifth condition, not required by the characteristic analysis, is 
used to close the set of equations. The need for this fifth condition could 
be removed by the use of one-sided differencing or by applying one of the 
governing equations at the wall. In this study for convenience the second 
method was used and the boundary layer approximation to the normal momentum 
equation was applied at the wall. This can be expressed as 


n w *V p =° (22) 

Studies nave indicated that there is little difference between using this 
equation and the nonapproximated normal momentum equation. 

The symmetry conditions are meant to be applied on a plane or axis of 
symmetry. The velocity conditions used require that the cross plane velocity 
component normal to the axis or plane of symmetry equals zero, i.e., 


n s • V = O 


(23) 


where n s is the unit vector normal to the axis or plane of symmetry and 
that the first derivatives of the remaining two velocity components equal 
zero. Two other conditions must be set on the axis or plane of symmetry. 
Usually the symmetry conditions on pressure and temperature are used, viz. 


n s -7P--0 


( 24 ) 


19 



and 


n s * VT = 0 (25) 

The final category of boundary conditions used in this investigation is the 
boundary conditions to be used on what can be considered to be external 
surfaces. Specifically these boundary conditions are meant to be used on the 
external boundary of the coflowing ambient stream. In this case a set of 
boundary conditions that will allow all disturbances which originate from 
within the computational domain to pass through the boundary without spurious 
reflection. The technique used in this investigation is based on the concept 
that in a simple wave region the flow properties remain constant along Mach 
lines. Thus the first derivatives of the flow variables in the direction of 
the Mach angle should be small and are here set equal to zero. The technique 
is termed Mach wave extrapolation and yields the boundary conditions 

rT m - V\7 = 0 (26) 

TT m VP=0 (27) 


and 




VT = 0 


(28) 


■> , , 
where n m is the unit vector in the direction of the local Mach angle. This 

technique requires computation of the Mach angle and has been successfully 

applied to a number of test cases both by the present authors and the authors 

of Ref. 25 and 26. 


I 
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SOLUTION OF THE GOVERNING EQUATIONS 


Numerical Methods 

The governing equations in both the supersonic and the embedded subsonic 
portion of the flow are simultaneously solved by the consistently split 
linearized block implicit (LBI) technique described in detail in 
Refs. 27 and 28. This technique can be logically divided into three parts: 
(1) linearization of the governing equations (2) discretization of the 
resultig set of linearized equation by finite difference approximation of 
derivative terras and (3) simultaneous solution of the resultant set of linear 
coupled algebraic equations. Application of the LBI technique to a set of 
governing equations (and boundary conditions) that is well posed for forward 
marching is straightforward. The system of governing equations can be 
written in the following form: 

<3H (<f>) (29) 

— = D(<£) + 


where 4> is the column vector of dependent variables (wj, W2, W3, p, h Q ) , 

H and S are column vector algebraic functions of $, and D is a column vector 
whose elements are the spatial differential operators. 

When a solution at the i + l sc station, at some distance AX 
downstream, is desired after a solution is obtained at an arbitrary i ^ 
streamwise station, the solution procedure is based on the following implicit 
marching direction difference approximation of Eq. (29) 


i + I i A i + I i + 

(H - H )/Ax = (D + S 


(30) 


where, for example, H 1+ ^ denotes H(<|> 1+ 1). A local spatial linearization 
(Taylor series expansion about <fr r ) of requisite formal accuracy is 
introduced, and this serves to define a linear differential operator L such 
that 

d' + 1 = D 1 + L* (<£' + l - <f>' ) + 0(Ax 2 ) (31) 
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Similarly 


h' + 1 . h' + (-§■)' 0(A«*) 

i + I i / d S \ i , , i + I . i , 2, 

S = s + ( _ ^ _ ) (< £ "0 > + 0 ( Ax 2 ) 


(32) 

(33) 


Eqs. (31-33) are inserted into Eq. (30) to obtain the following system which 
is linear in 4>i + l 


(A — AxL')(<£ i+l - 4 >') = Ax(0'+ S 1 ) (34) 

and which is termed the linearized block implicit (LBI) scheme. Here A 
denotes a square matrix defined by 



Eq. (35) is 0(AX) accuracy. 

To obtain an efficient algorithm, the linearized system, Eq. (34) is 
split using ADI techniques. To obtain the split scheme, the multidimensional 
operator, L, is rewritten as the sum of two 'one-dimensional* sub-operator 
L^(i = 2,3) each of which contains all terms having derivatives with 
respect to the i^-cross plane coordinate. The split form of Eq. (34) can 
be derived either as in Ref. 27 by following the procedure described by 
Douglas and Gunn (Ref* 29) in their generalization and unification of scalar 
ADI schemes, or using the approximate factorization as in Ref. 30. For the 
present system of equations, the split algorithm is given by 

(A - AxL|) (<£*-<£') = Ax(d‘ + S 1 ) (36) 

(A- AxLg) (<£' +l -<£') = (37) 

where <J>* is the consistent intermediate solution (Ref. 28). If spatial 
derivatives In and D are replaced by the difference formulae indicated 
previously, then each step in Eq. (36) and Eq. (37) can be solved by a block 
tridiagonal elimination. 
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Combining Eqs. (36) and (37) gives 


(A-AxLj)A"^A-AxLj,)(<£ i + l -<£ i ) = Ax(D‘ + s') 

which approximates the unsplit scheme, Eq. (34) to 0(AX ). Since the 
intermediate step is also a consistent approximation to Eq. (34), physical 
boundary conditions can be used for <f>* (Refs. 28, 31). 


I 

! 
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TEST CASES 


A primary objective of this study was to develop a computational 
procedure which is capable of analyzing a supersonic jet flow interacting 
with an external stream. Therefore, three three-dimensional and one 
axisymmetric test cases were considered to demonstrate the capability of the 
analysis and the associated computer program (PEPSIN). Three-dimensional 
cases deal with the underexpanded supersonic jet flow (i.e., Pj > Pe) 
exhausting from the rectangular nozzle of three different aspect ratios (1, 2 
and 5 respectively). As discussed in the Introduction Section, a supersonic 
jet was assumed to interact with a supersonic external stream. However, 
since no data were available to validate the results obtained by PEPSIN 
analysis for these three-dimensional cases, it was decided to consider an 
axisymmetric case (Figure 3) first because of the existing computational data 
(Ref. 17) to validate the PEPSIN results. Compared to the three-dimensional 
case, the axisymmetric results are, in general, easy to interpret. Moreover, 
experience and insight obtained during the study of the axisymmetr ic case is 
useful for the study of three-dimensional cases. Most of the flow conditions 
utilized in the computation of the axisymmetr ic case were obtained from 
Ref. 17. Some of the unavailable conditions were either assumed or 
estimated. The analysis was performed in the same axisymmetric geometry as 
reported in Ref. 17. The Mach number of both jet and external stream was 
2.0. The jet flow was hot (T j = 1500°K) and the static pressure ratio of 
jet with respect to the external stream was Pj/P£ = 1.45. The static 
temperature of the external stream was 240°K. The static pressure of the 
external flow was assumed as 2864 (Newton/m ). Thus, the computed Reynolds 
number per unit length based on the ambient flow conditions was 
1.6685 x 10 6 . To analyze the interaction of the jet flow with the coflowing 
ambient supersonic flow, 98 mesh points were distributed in the radial 
direction. The grid points were tightly packed in the vicinity of 
estimated turbulent mixing layer region to achieve the maximum resolution. 

The continuity equation, two momenta equations and the energy equation were 
used as governing equations. For the present case, the subsonic regions 
associated with the nozzle boundary layer and the base of nozzle wall are 
embedded in the supersonic flow. However, these embedded subsonic regions 
appeared to be excluded from consideration in Ref. 17 because the technique 
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used in Ref. 17 was not capable of analyzing such regions. Therefore, since 
it was desired to exercise the PEPSIN code under the same conditions as in 
Ref. 17 for this test case, the solution was marched with the initial 
profiles devoid of subsonic flow in this study. However, it should be noted 
that proper treatment of the subsonic regions may be very important in 
simulating physical situations. Since Ref. 17 didn't specify explicitly how 
the initial profiles were generated, the shape of initial profiles in the 
present calculations was approximately represented by the profiles which 
would be expected downstream in the mixing layer. For instance, the 
resulting velocity profile at the initial station would consist of two 
different uniform streams connected by the smooth hyperbolic tangent curve 
which approximate the transverse distribution of streamwise velocity 
component in the mixing layer. The boundary conditions imposed on the axis 
of symmetry are symmetry conditions, while the conditions on the external 
boundary are Mach line extrapolations. The streamwise step size used for 
marching the solution was 0.005 of the transverse size of computational 
domain, i.e., 0.025 of jet radius (corresponding to a Courant number of 
1.26). 

The computed results by the PEPSIN code are compared with the 
calculation of Dash et al * (Ref. 17) in Figures 4(a) and 4(b). The figures 
show static pressure variation in the streamwise direction at two different 
radial locations (along the jet axis and r/rj = .5 where rj refers to the 
radius of axisymmetric jet). Considering the uncertainties of the initial 
flow conditions which were used in PEPSIN analysis, the agreement between two 
different calculations is reasonable. The PEPSIN analysis indicates a 
stronger streamwise damping of wave strengths compared to the computation 
(Ref. 17). This appears to be associated with the difference in turbulent 
dissipation which occurs as a result of using different turbulence models in 
the two separate computations. The PEPSIN analysis is based on a algebraic 
mixing-length turbulence model, whereas a two-equations model ( <e model) was 
utilized by Dash and co-workers. However, both calculations show the absence 
of any noticeable wave structure beyond the first shock cell, which indicates 
that the reflected wave resulting from the interaction at the end of the 
first shock cell becomes very weak. The favorable comparison with the 
results of Ref. 17 indicates that both PEPSIN approach and its computer 
program are fundamentally sound. 
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After the PEPSIN code was successfully validated for an axisymmetric 
case, three-dimensional cases were considered. As discussed above, the 
underexpanded supersonic jet was assumed to be exhausted from a rectangular 
nozzle of aspect ratio 1,2 and 5 respectively (Figure 5). An aspect ratio of 
a cross-section is defined by the ratio of spanwise dimension divided by the 
transverse dimension. Unlike the previously considered axisymmetric case, 
estimation of initial profiles for three-dimensional cases- was not 
straightforward. Based upon our experience with the previous supersonic flow 
calculations, an initial profile must be not only smooth but also consistent 
with the governing equations. If not, perturbations, inconsistencies, etc. 
can persist far downstream. The difficulty of being consistent with the 
governing equations has been found to be further pronounced for 
three-dimensional cases. Therefore, in this study the problem was slightly 
modified such that initial profiles could be generated more easily. Thus, 
the jet was assumed to be exhausted from a short straight duct with a square 
or rectangular cross-section. In the present study we are primarily 
concerned with developing a computer code capable of predicting the 
interaction of a supersonic jet with an ambient supersonic flow downstream of 
exit plane. Since no specific nozzle geometry was defined in the present 
investigation, the length of duct was kept short mainly to reduce the 
computation time. The PEPSIN code is, however, capable of analyzing the flow 
field within a longer nozzle of more complex geometry as shown by the 
applications of the PEPSIS code to inlet flow systems (Ref. 16). 

In this study the streamwise marching solution procedure began at a 
location upstream of the duct entrance. As an initial profile for this 
problem, the flow was assumed to consist of two coflowing uniform streams at 
different velocities as shown in Figure 5. One stream at faster speed moves 
into the duct, whereas the external flow moving at a slower speed moves 
around the duct. The initial conditions such as velocity ratio of two 
streams and corresponding pressure ratio were assumed to be the same as the 
previously considered axisymmetric case. As a consequence of the 
modification, the flow at the duct exit was expected to contain some extra 
physical features due to the development of flow near the duct entrance. The 
flow both within and outside the duct was complicated because of the shock 
waves generated due to the initial shock wave-boundary layer development on 
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the walls of the duct. In particular, the pressure field within the 
cross-section became non-uniform as one would expect in a rectangular duct 
flow. 

The Mach number of both the initial streams was 2.0, while the Reynolds 
number per unit length based on the free stream (external flow) properties 
was 1.6685 * 10^ (a jet exit Reynolds number was 2.847 x 10^). In order 
to analyze both the internal and external flow of the duct as well as the jet 
flow interacting with the external stream, 50 x 50 mesh points for aspect 
ratio 1 and 50 x 80 mesh points for both aspect ratio 2 and 5 were 
distributed in each cross-sectional plane. Mesh points were packed in the 
region of the interaction. Because of the geometrical shape of both duct and 
jet flow, two planes of symmetry can be associated with this problem, and 
hence the flow need only be solved in the quarter plane. On the planes of 
symmetry, symmetry boundary conditions were utilized. On the other hand, 
three no-slip conditions for the velocity components, normal momentum 
equation and adiabatic wall conditions were imposed on the solid walls. On 
the external boundary of flow symmetry boundary conditions were used again. 
For the three-dimensional cases, the continuity equation, three momenta 
equations and the energy equation were used as governing equations. 

The turbulence model utilized in this study was an algebraic mixing length 
model. This model was selected primarily because it is convenient to use at 
this point (demonstration of capability of new computer code was one of the 
primary goals in this study) and, furthermore, it can save some computer 
resources (time and storage). However, if a demonstrated need for advanced 
turbulence model exists, this can be easily accommodated within the existing 
code framework. The mixing length was assumed to be proportional to the half 
width of mixing region and, therefore, constant at each streamwise location. 
In the present study, for the purpose of the turbulence model the spreading 
angle of the mixing region was assumed to be 3°. 

The computed results of three-dimensional cases are presented in Figs. 8 
through 16. Since experimental data is not available for direct comparison, 
discussion will primarily focus on qualitative description of the essential 
physical features predicted by the computation. To assist the reader in 
interpreting the plots of various variables, both Figs. 6 and 7 are 
presented. Figure 6 illustrates a sideview of the computational domain 
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on a vertical plane when the computational domain is cut through the duct 
vertically. Figure 7 provides a cross-sectional view of the domain. 

As discussed before, three-dimensional test cases were run for three 
different aspect ratios of rectangular cross-sec t ion . Therefore, discussion 
will first concentrate on the description of the results for an aspect ratio 
equal to 1. Then, the effects of larger aspect ratios on the flow field will 
be considered. 

Calculated results in the form of stagnation pressure and 
cross-sectional secondary velocity vectors for square nozzle case 
(aspect ratio 1) are presented in Figs. 8-9. Stagnation pressure isobars in 
the cross-sections at seven streamwise locations are presented to define a 
streamwise development of plume boundary (Fig* 8). Initial upstream contours 
were obtained at a distance of X/H =0.33 from the exit plane of the duct 
where H is a height of the duct. H is chosen to be 0.4. Meanwhile, final 
downstream contours were obtained at a distance of X/H = 22 from the exit 
plane of the duct. The plume boundary experiences a severely distorted 
square shape before evolving into a smooth shape far downstream. This is 
probably due to strong three-dimensional effects of the flow. Compared to 
the usual supersonic jet exhausted from a nozzle, these three-dimensional 
effects may have been more pronounced by the use of a short straight duct as 
a means of avoiding the difficulty associated with generating 
three-dimensional initial profiles. As discussed before, a uniform flow was 
assumed to move into the duct entrance. As a consequence, the flow within 
the duct becomes highly three-dimensional due to the generation of shock 
waves associated with the leading edge. Combined with the three-dimensional 
effects, a highly non-uniform pressure develops in the c ross-sec t ion of 
duct. In particular the pressure in the corner region is large. Thus, the 
largest pressure difference across the duct wall occurs in the corner 
region. Because of this non-uniform pressure distribution , a distortion of 
the plume boundary occurs. As can be seen in Fig. 8, the distortion is 
slightly asymmetric about the bisecting line (45° line) from the corner. The 
asymmetry was caused by the approximations used in the subsonic layer of the 
corner region of a nozzle. As discussed in the analysis section, both the 
continuity and normal momentum equations are approximated. The nearest solid 
wall from a grid point is identified to implement the subsonic 
approximations. However, the grid points located on the bisecting line have 
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two such surfaces. In this case the upper surface was arbitrarily chosen 
and, hence, a slight asymmetry introduced. It should be noted that 
asymmetry can be reduced by providing more grid points in the corner region, 
and as the flow proceeds downstream this effect significantly diminishes. A 
streamwise development of secondary velocity vectors in the cross-section is 
presented in Fig. 9. A strong secondary flow is seen in the corner region. 

As a result of the non-uniform pressure distribution, the flow develops two 
streamwise vortices which will eventually be diffused by the turbulent mixing 
action. To show the effects of jet aspect ratio on the flow field, Figs. 10 
and 11 are presented. Both figures illustrate a streamwise development of 
stagnation pressure contours. Initial upstream contours are obtained at a 
distance of X/H - 0.675 from the exit plane of the duct where H is a height 
of duct. H is taken to be 0.2. Final downstream contours are obtained at a 
distance of X/H 88 44. 

The outer edge of the mixing region may be defined qualitatively based 
on these two figures. Examination of the plots indicates that the mixing 
layer spreads out more quickly in the corner region as discussed previously 
in the case of aspect ratio 1. Streamwise development of cross-sectional 
plume boundary is very similar to that for aspect ratio 1, although 
distortion of plume boundary is localised in the corner region at an aspect 
ratio 5. It appears that the rapid mixing in the corner region is related to 
the three-dimensional effects. However, as both figures show, the plume 
boundary gradually approaches a smooth elliptic or circular shape in far 
field. In the near field downstream of the nozzle exit, the high pressure in 
the corner region associated with the shock waves may have played a 
significant role in enhancing the mixing. It should be noted that similar 
mixing behavior may not be observed when the calculation is performed for a 
realistic nozzle under the initial profiles based on the exit flow 
conditions . 
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DISCUSSION AND CONCLUSIONS 


The primary objective of this investigation was to develop and validate 
a numerical procedure for the calculation of the interaction of 
three-dimensional supersonic jet with an external flow. The new computer 
code (named PEPSIN) was developed by extending the capability of existing 
code (PEPSIS) suitable for an analysis of jet flows. The major difference 
between these two codes lies in the new capability of calculating flow fields 
when the computational domain contains internally embedded solid bodies. 

To validate the code, one two-dimensional and three three-dimensional cases 
were computed. Comparison of the two-dimensional results with other 
computations was performed to validate the new code for this class of cases. 
Three-dimensional calculations were successfully performed, and essential 
physical features were predicted. 
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USERS' MANUAL 


The PEPSIN Users' Manual is meant to serve as a guide in helping the 
user make successful runs with the PEPSIN computer program. The degree of 
success obtained by the user will depend on the skill of the user and his 
ability to correctly apply the code to his particular problem. The code will 
solve the governing equations, subject to the user supplied boundary 
conditions, however, meaningful results will only be obtained if the boundary 
conditions are appropriate to the problem. In addition, the user must 
specify viscosity models, initial conditions, a coordinate system and the 
location of grid points to adequately resolve the flow. The user with a good 
knowledge of the physics involved in his problem and how the code models the 
physics should with a moderate amount of experience be able to successfully 
apply the code to a wide variety of supersonic flow problems. 

The Users' Manual is divided into eight parts consisting of: (1) a flow 

diagram, (2) a brief description of each subroutine and its use, (3) a list 
of the Fortran variables and a description of their meaning, (4) a 
description of the logical file units utilized by the PEPSIN computer code, 

(5) a detailed description of the input required by the PEPSIN computer code, 

(6) a description of the common error conditions that may be encountered 
during the execution of a PEPSIN run and the corrective action to be taken, 

(7) sample input for a three-dimensional case and (8) sample output for the 
corresponding case. 


Flow Diagram 

The purpose of the flow diagram is to help the user understand the basic 
flow of information within the PEPSIN computer code. Because of the size of 
the code (approximately 13,000 cards), a detailed flow diagram would be 
prohibitively large and probably be of little value to the user. Therefore, 
the flow diagram is intended only to give a general overview of the structure 
of the code. The interested user is urged to consult the program listing for 
details . 


31 




32 


































PEPSIN Subroutines 


Subroutine 

ADDRES 

ADI 

ADICP 

ADIUN 

AMARCH 

AMATRX 

ARTVIS 

AVRG 

BC 

BLKDATA 

BLT 

BULEEV 

CONVCT 

CORBND 

CORTRN 

CROSEC 

CURVT 

DATAS 

DELTX 

DELTXZ 

DIFF 


Purpose 

Calculate addresses for finite difference 
representation of metric and fluid dynamic 
variables. 

Master control subroutine for ADL procedure. 

Control subroutine for coupled equations. 

Control subroutine for uncoupled equations. 

Linearizes streamwise convective terms. 

Linearizes all streamwise terms. 

Artificial dissipatior subroutine. 

Calculates averaged quantities in cross plane. 

Boundary condition subroutine. 

Stores default values of key variables. 

Calculates boundary layer thickness. 

Calculates Buleev mixing length. 

Linearizes cross plane convective terms. 

Calculates geometry transformation information on 
boundaries . 

Calculates geometry transformation information for 
interior points. 

Control subroutine for calculation of derived 
variables . 

Linearizes curvature terms. 

Logical file control subroutine. 

Calculates transformation information for ICORD = 2 
option. 

Calculates transformation information for ICORD = 3 
option. 

Linearizes diffusion terms. 


Subroutine 

Purpose 

DISFCN 

Calculates dissipation function. 

DIV 

Calculates divergence of velocity. 

D0P2 

Control subroutine for linearization of Y-direction 
and source terms. 

D0P3 

Control subroutine for linearization of Z-direction 
terms. 

ENDCAP 

Control subroutine for endcap conditions. 

EOS 

Equation of state subroutine linearizes and updates 
pressure and temperature. 

FGFUN 

Calculates geometry groupings. 

GAUSS 

Solves uncoupled tridiagonal set of equations. 

GENCBC 

Control subroutine for coupled boundary conditions. 

GENUBC 

Control subroutine for uncoupled boundary 
conditions. 

GEORD 

Controls setup of computational domain and 
calculates geometry at each cross-section. 

GEOTRB 

Generates ssetric information* 

INDIC 

Determines if flow is subsonic or supersonic at grid 
points. 

INPUTS 

Input subroutine. Input data enters and is 
processed. 

INTEBC 

Performs a two-dimensional linear interpolation for 
wall transpiration rates. 

LAMP 

Calculates laminar profile. 

LAW 

Calculates nondimens ional velocity, U + as a 
function of nondimens ional distance Y + . 

LENGTH 

Calculates mixing length. 

LOOP 

Determines the loop index at a grid point in both Y 
and Z-direction. 

MAIN 

Main control program. 

MATPRT 

Prints elements of block tridiagonal matrix. 
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Subroutine 


MGAUSS 

MGERR 

NMLIST 

OUTPUT 

PLOT 

PLOTIN 

PROF 

QUICK 

RE ADZ 

RESTRT 

ROTATE 

SETBVL 

SETBVP 

SETUP 

SHEAR 

SONIC 

SPREAD 

STORG 

SUB 

SWITCH 

TANHYP 


Purpose 

Control subroutine for solving block tridiagonal 
systems of equations. 

Calculates error associated with solving block 
tridiagonal system of equations. 

Subroutine for printing namelist input information. 

Control subroutine for printing out results on a 
cross-sectional (Y-Z) plane. 

Writes plot information on logical file unit JPLOT. 

Writes first record of general information on 
logical file unit JPLOT. 

Generates initial profiles. 

Matrix elimination subroutine. 

Prepares variables for printing. 

Reads and writes restart information. 

Rotates data from columns to rows and vice versa. 

Updates boundary information a line at a time. 

Updates boundary information a point at a time. 

Determines the extent of computational domain, 
number of loops and subsections at each 
cross-section 

Control subroutine for the calculation of wall shear 
velocity. 

Determines the extent of supersonic and subsonic 
region. 

Spreads two-dimensional data to three dimensions. 

Control subroutine of the storage of temperature and 
density on boundaries for viscosity calculation. 

Contain special subsonic logic. 

Calculates streamwise location for switch of 
boundary condition. 

Grid stretch subroutine. 
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Subroutine Purpose 


TNDER 

Calculates normal derivative of temperature at a 
wall. 

TRANS 

Transition model subroutine. 

TURB 

Turbulence model subroutine. 

TURBP 

Calculates turbulent profile based on theory of 
Maise-McDonald . 

VISCOS 

Constant and laminar viscosity subroutine. 

WALLFN 

Calculates wall shear velocity. 

WHERE 

Determines the surface number of four nearest 
boundaries for a given point. 

WRMATR 

Writes block tridiagonal dump information on logical 
file device NUNERR. 

YCALC 

Calculates Y and Z locations. 

ZERO 

Zeros out linearization arrays# 
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Logical File Units Utilized by PEPSIN Computer Code 


The PEPSIN computer code utilizes up to eleven (11) logical file units 
during the execution of a run stream. In many cases not all eleven units are 
used, and hence in these cases there is no need to define all eleven units. 
All references to a logical file unit in the PEPSIN computer code is 
accomplished through the use of a FORTRAN name rather than through a specific 
unit number. Thus, if the user desires to change a logical file unit number, 
this can be done through the input file. A list of the logical file units 
utilized by the PEPSIN computer code, their FORTRAN name, default value unit 
number, and a brief description of the use of the unit is presented below. 

All units are sequential. 


FORTRAN Name 


Default Unit Number 


Description 


MIN 

MOUT 

MASS1 


MASS2 


MSDD 


JDRUM 


5 Input data unit. 

6 Printed output unit. 

8 First unit which stores 

dependent and derived 
variables either by rows or 
columns. Not needed for 
two-dimensional cases, i.e., 
when TWOD = .TRUE. 

9 Second unit which stores 
dependent and derived 
variables either by rows or 
by columns. Not needed for 
two-dimensional cases, i.e., 
when TWOD = .TRUE. 

15 Unit which stores dependent 

and derived variables by rows 
or columns. Not needed for 
two-dimensional cases, i.e., 
when TWOD - .TRUE. 

11 Unit which contains output of 

ADD computer code. Only 
needed when IGEOM = 10 or 11. 
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FORTRAN Name 

Default Unit Number 

Description 

KDRUM 

12 

Unit which stores final 
metric information. Needed 
in all cases. 

NUNERR 

14 

Unit which stores information 
concerning the block 
tridiagonal matrix 
inversion. Needed when MGDMP 
* 0 . 

JPLOT 

15 

Unit which stores plotting 
information. Needed when 
IPLOT * 0. 

JRSTIN 

10 

Input restart unit. This 


unit contains appropriate 
common block information and 
the values of the dependent 
and derived variables at each 
cross-sectional grid point at 
the restart streamwise 
station. Needed only when 
IRSTIN * 0. 

JRSTOT 10 Output restart unit. This 

unit contains appropriate 
common block information, and 
the values of the dependent 
and derived variables at each 
cross-sectional grid point at 
the restart streamwise 
station. Needed only when 
IRSTOT * 0. Default is 
JRSTIN = IRSTOT; however, it 
is desired to have separate 
input restart and output 
restart files set JRSTOT = 

17. 
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PEPSIN Input 


Except for an initial title card and plot file input data, the entire 
PEPSIN input is entered by means of the NAMELIST format. There are two 
primary advantages to the use of the NAMELIST format: (1) that if the 

default values (defined in the block data subroutine) are acceptable, the 
user need not input that variable, and (2) the order (within a given 
NAMELIST) in which the variables are entered is irrelevant. There are seven 
NAMELIST input files in the PEPSIN code, $REST, $LIST1 through $LIST5 and 
$LISTR. The first file is read in the main program and enters restart 
information. The second through sixth NAMELIST files are read in subroutine 
INPUTS. Basically, the NAMELISTS $LIST1 through $LIST5 can be divided by 
function. $LIST1 enters information about the governing equations and 
appropriate boundary conditions, $LIST2 enters reference and free stream 
conditions, $LIST3 enters geometric information, $LIST4 enters viscosity 
model and initial profile information and $LIST5 enters file output 
information. The last NAMELIST file is read in subroutine CROSEC. 

$LISTR enters information needed to reset the grid array indicator, IFBW, 
for the computational domain in the cross-section. The PEPSIN computer code 
evolved from the PEPSIS computer code (Ref. 16) which was developed 
primarily for flows in supersonic inlets and for external flow situations. 

The primary difference between the two codes is in the ability of PEPSIN to 
consider cross-sectional geometric configuration that contain re-entrant 
corners as for example, might occur in the geometry illustrated in Fig. 12. 

In concept, the extension of the PEPSIS procedure to such geometric 
configurations is straightforward. However, to incorporate the re-entrant 
corner logic into the code in a general manner adds considerable complexity 
to the computer code. To minimize the effect of the complexities on the user 
while at the same time not detract from the generality that may be neces- 
sary at the present time or in the future, a very general procedure was 
incorporated into the PEPSIN code. This procedure utilizes the concept 
of a "grid array indicator 11 (FORTRAN variable IFBW) to type the grid 
points in the cross-sectional plane. Referring to Fig. 12, it can be 
seen that the idea is to type the grid points according to the function 
they serve. The convention is to type grid points as follows: 
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IFBW=1 

Indicates 

IFBW=2 

Indicates 


thus does 

IFBW=3 

Indicates 

IFBW=4 

Indicates 

IFBW-5 

Indicates 


that a grid point is an interior fluid point. 

that a grid point is an interior solid point, and 
not effect the calculation. 

that a grid point is a noncorner boundary point, 
that a grid point is an inward corner point, 
that a grid point is a re-entrant corner point. 


The FORTRAN variable is dimensioned IFBW(NN,NN). Based upon the information 
supplied by the user, subroutine SETUP sets up the necessary parameters to 
control the ADI procedure. Initially, the grid array indicator is introduced 
through NAMELIST LIST3. As the solution is marched downstream, the values of 
IFBW will remain unchanged. However, when it is desired to change the 
cross-sectional grid structure this can be accomplished by requesting that 
NAMELIST LISTR be read. This is determined by the input FORTRAN variable 
NBRKX which is read in NAMELIST LIST 3. 


In addition to the explanation of the grid aray indicator, it is 
convenient at this point to explain the surface number (KSURF) convention and 
direction (IDIR) convention utilized throughout the PEPSIN computer code. 
Several input variables, e.g. IBOUND( KSURF, IDIR) require a knowledge of this 
convention to correctly input data. The boundary surfaces are numbered 
according to a surface number relative to a cross-sectional computational 
direction. The convention is to allow values of IDIR=1 and 2 to correspond 
to the y and z physical directions, respectively. Referring to Fig. 12, it 
can be seen that the surfaces (with respect to y-direction) are numbered 1Y 
to 10Y. This corresponds to values of KSURF ranging from 1 to 10 and a value 
of IDIR corresponding to the y-direction i.e. IDIR=1 . The analogous 
convention is also applied to the z-direction surfaces. Thus, the value of 
KSURF represents the surface number relative to a given cross-sectional 
direction, IDIR. 

A description of all the PEPSIN input information will be given below. 


Card 1 

Columns 

Format 

Variable 

Function 

1 - 24 

6A4 

TITLE(I) 

Title Card 

Card 2 

Columns 

Format 

Variable 

Function 

1 - 2 

112 

ISYM 

Reciprocal of Symmetry 

3-12 

1F10.0 

SYSTEM 

SYSTEM = 1 - Quasi- 
Cartesian Coordinates 

SYSTEM = 2 - Quasi- 
Cyl indr ical Coord inatei 
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Namelist Input Description 


Namelist or 
Variable Name 

REST 


IRSTIN 


IRSTOT 


JRSTIN 

JRSTOT 


NFILE 


NSAVED 


Description 
Restart Options 

Marching station number where data is to be read in 
for restart case. 

IRSTIN = 0: Dead start case. 

IRSTIN * 0: Restart case started at station IRSTIN. 

Default: 0. 

Interval for saving restart information. 

IRSTOT = 0: No restart information is saved. 

IRSTOT * 0: Information is saved at each IRSTOTth 

station. 

Default: 0 

Logical file name of input restart file. 

Default: 10. 

Logical file name of output restart file. JRSTOT 
and JRSTIN do not have to be same file. 

Default: 10: 

File number on unit .JRSTIN desired for restart. 
Default: 0. 

Number of restart stations saved on JRSTOT. 

On a restart by setting JRSTOT = JRSTIN and NFILE - 
NSAVED, one file can be used for both reading and 
writing without destorying the information 
previously saved. 

Default: -1. 
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Namelist Input Description 


Namelist or 
Variable Name 

REST 


ICOMP 


LIST1 


IHSTAG 


IBOUND ( KSURF , IDIR ) 


IEQBC (KSURF , IDIR, IEQ ) 


JEQBC (KSURF , IDIR, IEQ ) 


subscript 


Description 
Restart Options 

Flag for computer options: 

ICOMP = 1: Uni vac computer option. 

ICOMP = 2: CDC computer option. 

ICOMP = 3: IBM computer - virtual memory option. 

ICOMP = 4: Disk writing computer option. 

Default: 4. 

Equations and Boundary Conditions 

IHSTAG = 0: Energy equation formulated in terms of 

static enthalpy. 

IHSTAG = 1: Energy equation formulated in terms of 

stagnation enthalpy. 

IHSTAG = 2; Stagnation enthalpy is constant. 
Default: 1. 

Boundary characteristics (wall or non-wall) 
on surface KSURF in a computational domain. 

IBOUND (KSURF, IDIR) = 1: Solid wall boundary. 

IBOUND (KSURF, IDIR) = 2: Non-wall boundary. 

Default: 40*1. 


Boundary condition of the governing equation IEQ 
on solid surface KSURF, 

Default: 120*2, 40*16, 40*11, 80*2. 

Boundary condition of the governing equation IEQ 
on non-wall surface KSURF. 

Default: 40*11, 20*2, 40*11, 20*2, 160*1. 

Boundary condition options used for either IEQBC or 
JEQBC are as follows: 

<t>: Any dependent variable. 

P: Pressure. 

T: Temperature, 

n = Normal to boundary. 

n = Unit vector perpendicular to axis of symmetry. 

C = Cartesian component. 
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Namelist Input Description 


LIST1 


Equations and Boundary Conditions 


Index Description 

1 = 0 (No change of <f> at boundary) 

2 <f> = 0 


3 


V or W known 


4 

5 

6 

7 

8 

11 


12 

13 

14 

15 


16 

17 


18 


pV or pW known 
AP = 0 

P = PRESS (KSURF,IDIR) 

AT = 0 

T = TWALL ( KSURF , IDIR ) 

3«f> ■ 0 (gradient of <|> normal to 
3n boundary specified at 

boundary) 

Mach Line Extrapolation using 
one-sided difference weights 

Slip boundary condition for 
velocity using wall function 

3P ■ 0 (gradient of pressure normal 
3n to boundary) 

3P = curvature (pressure gradient 
3n normal to boundary with 
curvature effects) 

Momentum equation in direction 
normal to boundary 

3T = 0 (adiabatic condition for 
3n wall or symmetry condition 

for non-wall) 

3T = DTDN ( KSURF, IDIR) 

3n 


19 Wall function boundary condition 
for temperature 

20 = 0 (same as 11, but applied at 

3n one grid point off the 
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Namelist Input Description 


LIST1 


Equations and Boundary Conditions 


Index 


Description 


21 


41 

42 

43 

44 

45 

46 

47 


Mach Line Extrapolation using 
central difference scheme at one 
point off the boundary 

8 2 4 > = 0 

an 2 

9 2 P = 0 
3n 2 

a = o 

an 2 

au„ = o 

3n 

av^ = o 

3n 

n c oVP = 0 
n c oVT = 0 


TWALL(KSURF.IDIR) 


PRESS (KSURF,IDIR) 


DTDN( KSURF, IDIR) 


ASW( KSURF, IDIR) 
BSW(KSURF, IDIR) 
CSW( KSURF, IDIR) 
DSW(KSURF,IDIR) 


Specified temperature on surface KSURF. 

Default: 40*1.0. 

Specified pressure on surface KSURF. 

Default: 40*0.0. 

Specified temperature gradient on surface KSURF 
Default: 40*0.0. 

Coefficient of a cubic polynomial fit for KSURFth 
surface to determine the axial location where 
boundary characteristics on surface KSURF should be 
changed from wall to non-wall or vice versa, i.e., 
IB0UND( KSURF, IDIR) automatically changed. 
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LIST! 


i 


LIST2 


IUNITS 


LREF 

REPL 

MINF 

PINF 

PR 

PRT 

PZERO 


Namelist Input Description 

Equations and Boundary Conditions 

Description 

Default: 

ASW = 40*1. OE + 10 
BSW = 40*0.0 
CSW = 40*0.0 
DSW = 40*0.0 

Freest ream and Reference Conditions 

Sentinel for units. 

IUNITS - 1: English units 

IUNITS = 2: Metric units 

Default: 2. 

Reference length, (ft or meters) 

No default. 

Reynolds number per unit length (ft -1 or m -1 ) 

No default. 

Frees tream Mach number. 

No default. 

Freestream static pressure (lbf/ft or nt/m ). 

No default. 

Laminar Prandtl number. 

Default: 0.74. 

Turbulent Prandtl number. 

Default: 1.0. 

2 2 

Freestream stagnation pressure (lbf/ft or nt/m ) 
No default. 
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Namelist Input Description 


LIST3 


ICORD 


NE(IDIR) 


Geometric Options 

Flag for coordinate transformation. 

ICORD = 1: Conformal coordinates 

ICORD = 2: Nonorthogonal coordinates 

X -► X 

Y -*• Y 

Z + c(x,z) 

ICORD = 3: Nonorthogonal coordinates 

X -*■ X 

Y ->• n(X,Y.Z) 

Z Z 

Default: 1. 

Number of grid points in the Y(IDIR = 1) and 
Z ( IDIR » 2) directions. 

No default values . 


NS Last streamwise station in each run. 

(Solution is marched from IRST1N+1 to NS). 

No default. 


XENTR, DELX, 

IAP(ICOUNT), AP(ICOUNT), 
DXMIN(ICOUNT), 

DXMAX ( I COUNT ) 


XENTR is the initial streamwise location. DELX is 
the initial stepsize in the streamwise (marching) 
direction, i.e., X(2) = XENTR + DELX. At streamwise 
station I the streamwise position is given by 
X(I) = X(I-l) + AP(X(I-1 ) - X(I-2) ) where if AP is 
greater than 1.0, the streamwise step size will 
increase by (AP-1.0) percent each step. If AP is 
less than 1.0, the streamwise step size will 
decrease by (1.0-AP) percent each step. DXMIN and 
DXMAX are lower and upper overriding limits on the 
step size. AP, DXMIN and DXMAX are dimensional so 
that streamwise step size variation can be changed 
by the LAP parameter, IAP denoting the streamwise 
location where these variables change. Values of 
XENTR and DELX should normally only be set on the 
initial run as these variables are automatically 
calculated for restarts. Maximum ICOUNT is 10. 
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Namelist Input Description 


LIST3 


XENTR, DELX, 

IAP(ICOUNT, AP(ICOUNT), 
DXMIN ( ICOUNT ) , 
DXMAX(ICOUNT), 
(CONTINUED) 


IGEOM 


TWOD 


TTl(IDIR) 


TT2(IDIR) 


Geometric Options 


Default : 


XENTR = 0.0 
IAP = 1,9*1000000 
AP = 10*1.0 
DXMIN = 10*0.0 
DXMAX = 10*1. 0E + 06 

DELX = No default. 


Flag for coordinate options 

IGEOM - Is Cartesian coordiates 

IGEOM = 2: Cylindrical coordinates 

IGEOM « 3: Polar coordinates 

IGEOM = 10: General orthogonal coordinates 

(Cartesian In cross plane) 


IGEOM -11: 


uerauxt; ±. 


General orthogonal coordinates 
(axisymmetric) 


Sentinel for two-dimensional option. 

If TWOD = .TRUE. TWO DIMENSIONAL 

If TWOD = .FALSE. THREE DIMENSIONAL 

Default: .FALSE. 

Grid distribution factor (lower surface (IDIR=1) - 
or, left surface (IDIR*2)). The closer the value is 
to 1.0, the tighter the packing. 

Default: 2*0.0. 

Grid distribution factor (upper surface (IDIR=1) - 
or, right surface (IDIR=2)). The closer the value 
is to 1.0, the tighter the packing. 

Default: 2*0.0. 
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Namelist Input Description 


LIST3 


T2( IDIR) 


XO(IDIR) 


YS(ILIM,IDIR) 


IFBW(JY,KZ) 


XBRKX( I COUNT ) 


Geometric Options 

Grid distribution factor (Y-direction (IDIR=1) - or 
Z -direction (IDIR= S 2)). The larger the value is, the 
tighter the packing. Use in conjunction with XO. 

The value should not exceed 5.0. 

Default: 2*0.0 

Clustering location of grid points within the 
computational domain (X0(1) is Y-location and X0(2) 
is Z-location). Use in conjunction with T2. 

Default: 2*0.0 

Define computational domain in Y-Z cross plane. 

ILIM defines either lower or upper limit, 

ILIM-1 (LOWER), ILIM=2(UPPER) • 

YS(1,1) - 0.0 - lower limit Y-direction 
YS(2,1) * 1.0 - upper limit Y-direction 
YS(1,2) = 0.0 - lower limit Z-direction 
YS(2,2) ■ 1.0 - upper limit Z-direction 
Default: 0.0, 1.0, 0.0, 1.0. 

Grid array indicator needed to set up the 
computational domain in the cross plane at each 
s t reamwise s t at ion . 

IFBW( JY,KZ) Notation 

IFBW = 1 - Interior fluids point 

IFBW = 2 - Interior solids point 

IFBW 85 3 - Conventional boundary point 
(fluids or solids) 

IFBW - 4 - Boundary point outward corner 

IFBW = 5 - Boundary point inward corner 

No default. 

Streamwise (marching direction) physical locations 
where an embedded solid body starts or ends. 

Maximum ICOUNT is 5. 

Default : 5*1 . 0E+06 . 


52 



LIST3 


NBRKX 


IBRKMN ( KZ , ICOUNT ) 

i 

i 

! 

i 

1 IBRKMXOCZ, ICOUNT) 


I 

' IEDGE 


LIST4 


IBCP(LP) 


DELTAP(IBCP.LP) 


CFP(IBCP.LP) 


Namelist Input Description 


Geometric Options 


Total number of the streamwise locations where 
embedded solid bodies start or terminate. 

Default: 0. 

Minimum Y index of the grid points which change 
their characteristics at each spanwise (KZ) location 
of the cross-section. Maximum ICOUNT is NBRKX. 
Default: 500*0. 

Maximum Y index of the grid points which change 
their characteristics at each spanwise (KZ) location 
of the cross-section. Maximum ICOUNT is NBRKX. 
Default: 500*0. 

Flag which tells whether grid point changes its 
characteristics from fluids to solids or vice versa. 

IEDGE = 0: From fluids to solids. 

IEDGE = 1: From solids to fluids. 

No default. 

Initial Profile, Turbulence Information 

Basic surface for initial profile generation for 
LP th loop. 

Boundary Layer Profile in a loop 

at lower surface 1 
at upper surface 2 
at left surface 3 
at right surface 4 

Default: 1. 

Boundary layer thickness on surface IBCP of LPth 
loop needed to generate the initial profile 
referenced to each surface. 

No default. 

Skin friction coefficient on surface IBCP of LPth 
loop needed to generate the initial turbulent 
boundary layer profile. 

No default. 


IBCP = 1 
IBCP = 2 
IBCP = 3 
IBCP = 4 
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IPROF 


Flag for initial profile options. 


IPROF = 
IPROF = 
IPROF = 

IPROF = 
Default 


1: Freestream profiles 

2: Initial profiles supplied by user 

3: Boundary layer profiles based on 

necessary input 

4: Same as IPROF^S, but angular components 

are obtained for general orthogonal 
coordinates. 

value is 1 . 
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Namelist Input Description 


LIST4 Initial Profile, Turbulence Information 

IMIXL Flag for mixing length options. 

IMIXL = 1: McDonald-Camarrata mixing length model 

based on prescribed boundary layer 
thickness (DELTAB). With wall shear 
value used to calculate nondimens ional 
distance 

j IMIXL = 2: Buleev mixing length model 

I IMIXL = 3: McDonald-Camarrata mixing length model 

i based on dynamically obtained boundary 

! layer thickness with fixed wall shear 

| IMIXL = 4: Same as IMIXL * 1, but local shear is 

used to calculate nondlmensional 
distance 

| IMIXL ■ 5: Same as IMIXL ■ 3, but local shear is 

used 

IMIXL = 6: Mixing length model for jet flow 

Default: 1. 

' BETA Angle of attack in degrees. 

Default: 0.0. 

YAW Yaw angle in degrees. 

Default: 0.0. 

IVISC Flag for viscosity options. 

IVISC * 1: Constant viscosity 

IVISC = 2: Laminar viscosity obtained from 

Sutherland's relation 

IVISC = 3: Turbulent viscosity is obtained 

from mixing length model 

IVISC = 4: Turbulent viscosity obtained from 

TKE - mixing length model 

Default: 1. 

DELTAB (KSURF, IDIR) Specified boundary layer thickness on surface KSURF 

for mixing length model of turbulence. 

No default values. 
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LIST4 


ITRANS 


IBLT 


TKEINF 

NPROF(LP) 

IDIRP 


LIST5 


IVARPR(I) 


Namelist Input Description 

Initial Profile, Turbulence Information 

Flag which tells whether transition turbulence model 
logic is used. 

ITRANS = 0: No transitional model is used 

ITRANS * 0: Transitional model is used 

Default: 0. 

Flag which tells whether boundary layer thickness is 
input or calculated dynamically. 

IBLT = 0: Boundary layer thickness is input 

IBLT t 0: Boundary layer thickness is dynamically 
calculated. 

Default: 0. 

Freestream turbulent kinetic energy. 

Default: 0.0. 

Number of initial profiles generated for each loop 
in Y or Z-direction. 

Default: 10*0. 

Basic direction for initial profile generation. 

IDIRP = 1: Y-direction 

IDIRP = 2: Z-direction 

Default: 2. 


File Output Information 

Index of variables to be printed. Needed only 
for three-dimensional flow. 

IVARPR(I) =0: No print 

IVARPR(I) = 1: Print every IPRINT steps 

IVARPR(I) = 2: Print every JPRINT steps 

I = 1: UVEL 

1=2: WEL 

1=3: WVEL 

I = 4: Density 

I = 5: Enthalpy 

I = 6: Turbulent kinetic energy 

I = 7: Turbulent dissipation 

Default: 5*1, 2*0, 3*1, 7*0, 1, 3*0, 1, 3*0. 
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Namelist Input Description 


LISTS 


IPLOT 


Initial Profile, Turbulence Information 


I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


8: Pressure 

9: Temperature 

10: Mach No. 

11: Mach No. indicator 

12: Stagnation temperature 

13: Stagnation pressure 

14: Pressure coefficient 

15: Laminar viscosity 

16: Mixing length 

17: Turbulent viscosity 

18: Effective viscosity 

19: Dissipation function 

20: Cell Reynolds number in Y-direction 

21: Cell Reynolds number in Z-direction 

22: ISS and JBOUND 

23: Heat transfer coefficient, skin friction 

coefficient and heat transfer rate 
24: Boundary layer thickness 

25: Cross-sectional average of flow properties 


Marching station interval for storage of plotting 
information. 


IPLOT * 0: No plotting 

IPLOT * 0: Store plotting information every 

IPLOT station 


Default: 0. 


IPRINT Primary marching station interval for printing. 

Default: 1. 

JPRINT Secondary marching station interval for printing. 

Default: 1. 
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Namelist Input Description 


Reset Information for Cross-Section 


LISTR 


IFBW 

See 

LIST3 

I BOUND 

See 

LIST1 

IEQBC 

See 

LIST1 

JEQBC 

See 

LISTl 

DTDN 

See 

LIST1 

TWALL 

See 

LISTl 

PRESS 

See 

LISTl 
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Error Conditions in the PEPSIN Computer Code 


Failure of the PEPSIN computer code to successfully execute a runstream 
can occur because of either inconsistent or incorrect input data or because 
of an attempt to apply the PEPSIN code to a case where the physics violate 
the assumptions inherent in the code. This section will address only the 
former mode of failure. Avoidance of the latter failure mode is dependent 
primarily on the users understanding of the basic physics of the case he is 
going to run, and the degree to which the PEPSIN code can be expected to 
model the physics. 

One method of discussing the inconsistent or incorrect input data mode 
of failure is by examining the possible failures in the various subroutines. 
Since the individual subroutines are responsible for separate tasks during 
the execution of a run, (e.g. overall control of the program geometry 
generation, etc.), this technique will in essence outline the possible 
failure modes as the tasks are performed. Discussion will occur in the same 
order as the run is executed. 

SUBROUTINE RESTRT 

There are two modes by which SUBROUTINE RESTRT can fail. Both involve 
improper use of the restart file. A message, RESTART INFORMATION REQUESTED 
AT (IRSTIN marching number) BUT STORED INFORMATION AT SEQUENCE (NFILE) IS AT 
STATION (Station number). This message occurs because the marching station 
number read off the NFILEth restart file does not match the input value of 
IRSTIN. The corrective action is to make NFILE and IRSTIN consistent with 
each other. Another possible mode of failure occurs when NFILE exceeds the 
number of files on the restart device, JRSTIN, in which case an END OF 
INFORMATION (or analogous statement) will appear in the day file. The 
corrective action is to recheck the input value of NFILE. If JRSTIN * 

JRSTOT, the value of NFILE is the number of the restart on device JRSTIN. 

SUBROUTINE INPUTS 

There are two failure modes in SUBROUTINE INPUTS. On the first case, 
the message NS = (input value of NS) GREATER THAN NSMAX = (dimension of X 
vector) will be printed if the number of marching stations exceeds the 
dimensional value of X, the streamwise locations. The corrective action is 
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to lower the value of NS. The second failure mode occurs when the Buleev 
turbulence model is specified for a two-dimensional case. Since this model 
is not applicable to two-dimensional cases, the message CANNOT USE BULEEV 
TURBULENCE MODEL IN TWO-DIMENSIONAL FLOW is printed. The corrective action 
is to specify an alternate turbulence model. 

SUBROUTINE SETUP 

If grid point characteristics indicator, IFBW, is not correctly 
specified, SUBROUTINE SETUP can fail. When the message FAILURE IN SETUP DUE 
TO INCORRECT IFBW is printed, the corrective action is to re check IFBW 
specified in INPUT DATA. IFBW should contain an information on the 
computational domain boundary which Is specified by boundary indicator. 

SUBROUTINE GEOTRB 

At present SUBROUTINE GEOTRB is coded to calculate metric information 
for values of IGEOM ■ 1, 2, 3, 10 and 11. Values of IGEOM 4-9 are left for 
various coordinates that may be coded In the future. Input value of IGEOM « 
4-9 will result in the message INVALID OPTION IN GEOTRB. The corrective 
action is to either change the value of IGEOM or to code In a new option. 

For IGEOM options 10 and 11 (conformal-Cartesian cross-section and 
conf ormal-axisymmetric cross-section) the metric information is externally 
generated by the ADD computer code. In this case, logical file units JDRUM 
and KDRUM must be defined. JDRUM contains the ADD code data which is then 
interpolated onto the PEPSIN mesh system. If the PEPSIN values of the 
streamwise coordinate is less than the first value of the ADD code streamwise 
coordinate no streamwise interpolation is possible and the message FAILURE IN 
GEOTRB - SQ12 = (PEPSIN position) SQ1 = (first ADD code position) SQ2 = 
(second ADD code position). The corrective action is to increase the value 
of XENTR (the first PEPSIN position) to a value greater than SQ1. On the 
other hand, if the value of a PEPSIN streamwise coordinate exceeds the last 
streamwise position generated by the ADD code and END OF INFORMATION message 
will appear In the day file. The corrective action is to either rerun the 
ADD code such that the maximum PEPSIN streamwise coordinate does not exceed 
the maximum ADD code streamwise coordinate or to reduce the maximum PEPSIN 
streamwise coordinate to an acceptable value. 
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SUBROUTINE INTEBC 


SUBROUTINE INTEBC performs a two-dimensional linear interpolation of 
the transpiration schedules on both X-Y planes at X-Z planes. If the number 
of st Teamwise stations on a surface at which data is input exceeds 15 (the 
dimensional size of the data arrays) a message FAILURE IN INTEBC VALUE OF 
NPTSX (surface number IBC) = (value of NPTSX(IBC) EXCEEDS DIMENSION LIMITS OF 
15 Is printed. The corrective action is either to updimension NPTSX and 
associated variables or to decrease the value of NPTSX. Likewise, in the Y 
or Z direction data can be input at up to 15 locations. If the value of 
NPTSYZ exceeds 15, the message FAILURE IN INTEBC VALUE OF NPTSYZ (streamwise 
location, surface number) = (value of NPTSYZ) EXCEEDS DIMENSION LIMITS OF 15 
Is printed. The corrective action is either to updimension NPTSYZ and 
associated variables or to decrease the value of NPTSYZ. 

SUBROUTINE QUICK 

If the choice of boundary conditions is incorrectly made, it Is 
possible that a singular matrix will result. This will manifest itself in 
SUBROUTINE QUICK in an attempt to divide by zero. The corrective action is 
to re-evaluate the choice of input boundary conditions to determine the 
source of the singularity. An example of an Improper choice of a boundary 
condition set would be to choose as boundary conditions the three no— slip 
conditions for the three momenta equations, the normal pressure condition for 
the continuity equation and the normal momentum equation for the enthalpy 
equation. In this case, the enthalpy does not appear in any of the boundary 
conditions, and hence a singular matrix would result. 

SUBROUTINE CROSEC 

Often, if a case is not going to successfully run, the code will cease 
operation in SUBROUTINE CROSEC. This will occur because of the existence of 
a negative temperature in which case the Mach number calculation will fail In 
SQRT. There can be many reasons for this failure mode. Usually, however, it 
can be related to inadequate numerical resolution of the physical processes 
that are occurring. For instance, a lack of transverse grid points might 
lead to large oscillations in the pressure or too large a streamwise step in 
the region where a wall inclination is rapidly changing might result in a 
temperature becoming negative. Sometimes it is difficult to know a_ priori 
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what grid resolution is necessary for a given problem. Usually, 
experimentation with two-dimensional cases can ironically provide some guide 
lines for three-dimensional cases. This in addition with the users 1 overall 
experience with the code and his understanding of the physical processes will 
usually provide the means of resolving the above problem. 

SUBROUTINE WHERE 

There are two modes by which SUBROUTINE WHERE can fail. For the first 
case, the message PROBLEM IN SUBROUTINE WHERE MSEGY = m JX ^ \ JY=n 2 KZ=n 3 
will be printed where ni, n2 and n3 refers to x, y and z location 
respectively if MSECY is not correctly specified, i.e., not consistent with 
the number of Y-perspective subsections. The corrective action is to 
increase the value of MSECY corresponding to the number of Y-perspective 
subsections. The second failure mode occurs when MSECZ is not consistent 
with the number of Z-perspective subsections. In this case, the message 
PROBLEM IN SUBROUTINE WHERE MSECZ = m JX=u \ JY=n£ KZ=n 3 is printed with nj, 
n2 and n3 referring to x,y and z location respectively. The corrective 
action is to increase the value of MSECZ corresponding to the number of 
Z-perspective subsections. 
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PEPSIN FORTRAN Variables 


FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


ACON 

LAUU 

CONSTANT IN ARGUMENT OF EXPONENTIAL FUNCTION 
FOR TRANSITIONAL MODEL 

AG(NN f 9,2) 

OPER 

DIFFERENCE WEIGHTS IN PHYSICAL COORDINATES 

AGEO 

GEOM 

COEFFICIENTS OF POLYNOMIAL FIT FOR BOUNDARY SHAPE 

AG1H(9) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AG IP 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AG2D(9) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AHP(5) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AH10(5,9) 

FGCOM 

TEMPORARY STORAGE ARRAY OF METRIC INFORMATION 

AIE(NN,7) 

PRFILE 

INITIAL PROFILE ARRAY 

AM ( NCPLD , 3ANCPLB+1 ) CCOM 

UIILITY MATRIX USED IN BLOCK MATRIX INVERSION 

AMACRT 

SUPER 

MACH NUMBER CRITERION USED IN LOCATING SONIC LINE 

AN(NEQS,NN) 

LIN 

STORAGE FOR LINEARIZATION COEFFICIENTS OF 
X - DERIVATIVES 

AP(10) 

GEOM 

AMPLIFICATION RATE OF MARCHING STEP SIZE 

APLUS 

LAUU 

CONSTANT IN ARGUMENT OF EXPONENTIAL 
FUNCTION FOR VAN DRIEST DAMPING FORMULA 

ASU(20,2) 

BOUND 

COEFFICIENTS OF POLYNOMIAL FIT FOR SWITCHING THE 
BOUNDARY SURFACE TYPE 

AVISC(2,NEQS) 

vise 

COEFFICIENT USED IN ARTIFICIAL DAMPING 

BETA 

REF 

ANGLE OF ATTACK 

BGEO 

GEOM 

COEFFICIENTS OF POLYNOMIAL FIT FOR BOUNDARY SHAPE 

BLOCK1 ( IADB3) 


STORAGE EUR ADD CODE METRIC INFORMATION 
EQUIVALENCED 10 C(l,l,l) 

BL0CK2( IADD3) 


STORAGE FOR ADD CODE METRIC INFORMATION 
EQUIVALENCED 10 C ( 1 , 1 ,NN/2+ 1 ) 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

BLIH(NN,5,4) 

LAWW 

DYNAMICALLY DETERMINED BOUNDARY LAYER THICKNESS 

BSW(20,2) 

BOUND 

COEFFICIENTS OF POLYNOMIAL FIT FOR SWITCHING THE 
BOUNDARY SURFACE TYPE 

BWD 

LIN 

CRANK-NICHOLSON FACTOR 

BUBI 

LIN 

INVERSE OF BWD 

C(NN f NCPLB r NN) 
CDUM(NDIM) 

CCOM 

BLOCK DIAGONAL MATRIX ELEMENTS 

TEMPORARY STORAGE ARRAY TO ROTATE DATA FROM COLUMNS 


TO ROWS AND VICE VERSA - EQUIVALENCED TO C(l f l,l) 
NDIM = MZVAR A HLEVEL A NN 


CFP( 4 , 10) PRFILB SKIN FRICTION COEFFICIENT 

CGEO GEOM COEFFICIENTS OF POLYNOMIAL FII FOR BOUNDARY SHAPE 

CMUINF VISC CONSTANT IN TURBULENT VISCOSITY MODEL 

CONGEO(llANN) GEOM COORDINATE TRANSFORMATION INFORMATION 

CONVDR UNITS CONVERSION FACTOR IN GOING FROM DEGREES TO RADIANS 

CONVRD UN IIS CONVERSION FACTOR IN GOING FROM RADIANS TO DEGREES 

C00R(NN,4) GEOM PHYSICAL COORDINATE INFORMATION WITH RESPECT 

TO ABSOLUIE ORIGIN AT N+1ST STREAMWISE LOCAIION 
COORN(NN r 4) GEOM PHYSICAL COORDINATE INFORMATION WITH RESPECT 

TO ABSOLUTE ORGIN AT NTH STREAMWISE STATION 
CP INF FREE FREE STREAM SPECIFIC HEAT 

CPREF REF REFERENCE PRESSURE COEFEICIENI 

CPREFI REF INVERSE OF REFERENCE PRESSURE COEFFICIENT 

CRITU OPER CRITICAL VELOCITY USED FOR FLARE APPROXIMATION 

CSOLN(NCPLB f NN) CCOM SOLUIION TO BLOCK TRI-DIAGONAL MATRIX INVERSION 

CSW<20,2) BOUND COEFFICIENTS OF POLYNOMIAL FIT FOR SWITCHING THE 

BOUNDARY SURFACE TYPE 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


CTWO 

REF 

NUMERICAL CONSIANI IN GOVERNING EQUATION 

CXI(NDIFM) 

ADDR 

STORAGE FOR FIRST DERIVATIVE DIFFERENCE UEIGHTS 

CXXKNDIFH) 

ADDR 

STORAGE FOR SECOND DERIVATIVE DIFFERENCE UEIGHTS 

Cl(NN) 

CCOM 

SUBDIAGONAL MATRIX ELEMENTS 

C1SUIH 

vise 

COEFFICIENT IN SUTHERLAND'S LAW OF LAMINAR 
VISCOSITY 

C2(NN) 

CCOM 

DIAGONAL MATRIX ELEMENTS 

C2SUTH 

vise 

COEFFICIENT IN SUTHERLAND'S LAW OF LAMINAR 
VISCOSITY 

C3(NN) 

CCOM 

SUPERDIAGONAL MATRIX ELEMENTS 

C4(NN) 

CCOM 

VECTOR ELEMENTS 

D 

VAR 

INDEX FOR DIVERGENCE 

DELHAX(NEQS) 

EQN 

MAXIM IUM VALUES OF IHE DELTAS 

BELTAB<20,2) 

LAUU 

SPECIFIED BOUNDARY LAYER THICKNESS FUR 
MIXING LENGTH MODEL OF TURBULENCE 

DELTAP(4 f 10) 

PRFILE 

BOUNDARY LAYER THICKNESS 

BELX 

6E0M 

STEP SIZE IN MARCHING DIRECTION 

DGEO 

GEOH 

COEFFICIENTS OF POLYNOMIAL FIT FOR BOUNDARY SHAPE 

DIF0P(6) 

OPER 

DIFFERENCE WEIGHTS IN COMPUTATIONAL COORDINATES 

DL 

VAR 

STORAGE LEVEL OF DIVERGENCE IN 
GENERAL PURPOSE STORAGE 

DM1 (NCPLB, NCPLD) 

CCOM 

GENERAL PURPOSE STORAGE ARRAYS FOR BLOCK 
TRI-DIAGONAL HATRIX INVERSION 

DM2( NCPLB) 

CCOM 

GENERAL PURPOSE STORAGE ARRAYS FOR BLOCK 
TRI-DIAGONAL HATRIX INVERSION 

DM3(NCPLD,NCPLD> 

CCOM 

GENERAL PURPOSE STORAGE ARRAYS FOR BLOCK 


r 
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FORTRAN 

COMMON 

SYMBOL 

BLOCK 

DS 

VAR 

DSL 

VAR 

BSU(20,2) 

BOUND 

DTDN(20,2) 

BOUND 

DTDNW(NN,5,4) 

BOUND 

DX 

OPER 

DXI 

OPER 

DXMAX(IO) 

GEOM 

DXMIN(IO) 

GEOM 

DI (NEQS,NDIFM, NN) 

LIN 

D1L(NEQS, NDIFM) 

LIN 

D2(NEQS,NDIFM,NN) 

LIN 

E(NN,NCPLD, NN) 


EPS 

VAR 

EPSMWF 

LAUU 

F< 1 4 , NN ) 

ZPLOT 


DESCRIPTION 

IRI-DIAGONAL MATRIX INVERSION 
INDEX FOR DISSIPATION FUNCTION 
STORAGE LEVEL OF DISSIPATION FUNCIION IN 
GENERAL PURPOSE STORAGE 

COEFFICIENTS OF POLYNOMIAL FII FOR SWITCHING THE 
BOUNDARY SURFACE TYPE 

SPECIFIED TEMPERATURE GRADIENT NORMAL TO 
BOUNDARY 

STORAGE ARRAY FOR TEMPERATURE GRADIENT NORMAL 
TO BOUNDARY 

STEP SIZE IN X-DIRECTION 
INVERSE OF DX 

MAXIMUM MARCHING STEP SIZE 

MINIMUM MARCHING STEP SIZE 

STORAGE FOR LINEARIZATION COEFFICIENTS OE 

Y - DERIVATIVES 

STORAGE OF LINEARIZATION COEFFICIENT 
FOR TEMPERATURE AND PRESSURE COMPUTATION 
STORAGE FOR LINEARIZATION COEFFICIENTS OF 
Z - DERIVATIVES 

ARRAY USED IN MGAUSS ERROR CHECK - EQUIVALENCED 
TO C(l,l,l) 

INDEX FOR DISSIPATION OF TURBULENCE 
KINETIC ENERGY 

CONVERGENCE CRITERION FOR WALL FUNCTION FORMULATION 
TEMPORARY STORAGE FOR PLOT INFORMATION 
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FORTRAN 

SYMBOL 

FACLH(20,2) 

FG(8,3 f NN) 
GAM(3) 

GAMMA 

GC 

H 

HFORM 

HINF 

HREF 

HREFI 

IA 

IADDO(NBIFM) 

IABDP(NDIFM) 

IADDSl(NN) 

IADDS2(NDIFM,NN) 

IADDS3(NDIFM f NN) 

IADDS4 (NDIFMAA2, 

IADDS5(NDIFMAA2, 

IADD1 

IADD2 

IADD3 

IAD I 

IADIM1 

IAP(IO) 


COMMON 

BLOCK 

LAUU 

METRIC 

GEOM 

REF 
UNITS 
VAR 
REF 
FREE 
REF 
REF 
CCOM 
ADDR 
ADDR 
ADBRF 
ADDRF 
ADDRF 
NN) ADBRF 
NN) ADDRF 
ADD 
ADD 
ADD 
SHEEP 
SWEEP 
GEOM 


DESCRIPTION 

MULTIPLICATION FACTOR TO BOUNDARY LAYER 
THICKNESS 

STORAGE FOR HEIRIC COEFFICIENTS AND DERIVATIVES 

COORDINATE TRANSFORMATION COEFFICIENT FOR 

BOUNDARY CONDITIONS 

RATIO OF SPECIFIC HEATS 

GRAVITY CONSTANI 

INDEX FOR ENTHALPY 

HEAT OF FORMATION 

FREE STREAM ENTHALPY 

REFERENCE ENTHALPY 

INVERSE OF REFERENCE ENTHALPY 

INDEX REFERRING TO SUBDIAGONAL MATRIX ELEMENTS 

ADDRESSES IN THE OPPOSITE DIRECTION 

ADDRESSES IN THE PRIMARY DIRECTION 

ADDRESS FOR POINT LOGIC OF FLUID VARIABLES 

ADDRESS FOR Y DERIVATIVE OF FLUID VARIABLES 

ADDRESS FOR Z DERIVATIVE OF FLUID VARIABLES 

ADDRESS FOR MIXED DERIVATIVE OF FLUID VARIABLES 

ADDRESS FOR MIXED DERIVATIVE OF FLUID VARIABLES 

NO. OF GEOMETRIC VARIABLES USED IN ADD CODE 

NO. OF TRANSVERSE GRIDPOINT USED IN ADD CODE 

RECORD SIZE USED IN ADD CODE 

ADI SWEEP DIRECTION 

IAD I - 1 

MARCHING STEP INDEX AT WHICH AP , DXM IN , DXMAX 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 

ARE REINITIALIZED 

IB 

CCOM 

INDEX REFERRING TO DIAGONAL MATRIX ELEMENTS 

IBC 

BOUND 

INDEX FOR BOUNDARY SUREACE IDENTIFICATION 

IBCP(IO) 

PRFILE 

BASIC SURFACE FOR INITIAL PROFILE GENERATION 

IBCPLP 

PREILE 

BOUNDARY SURFACE NUMBER IN EACH LOOP WHERE 
INITIAL BOUNDARY LAYER PROFILE IS REQUIRED. 
IBCPLP = 1 OR 2 IF IDIRP = 2 
IBCPLP = 3 OR 4 IF IDIRP = 3 

IBLT 

LAWW 

FLAG WHICH IELLS WHETHER BOUNDARY LAYER 
THICKNESS IS INPUT OR CALCULATED DYNAMICALLY 

IBOUND(20,2) 

BOUND 

BOUNDARY SURFACE TYPE INDICATOR 

IBRKMX(NN f 5) 

BOUND 

UPPER BOUNDARY OF AN EMBEDDED SOLID BODY 
EXPRESSED IN TERMS OF GRID POINT INDEX ON EACH 
Y-COORDINATE LINE AI UP TO 5 STREAMWISE LOCATIONS 

IEDGE 

BOUND 

SENTINEL WHICH SPECIFIES EITHER LEADING EDGE 
WHERE SOLID BODY STARTS OR TRAILING EDGE WHERE 
SOLID BODY TERMINATES. 

IBRKMN(NN,5) 

BOUND 

LOWER BOUNDARY OF AN EMBEDDED SOLID BODY 
EXPRESSED IN TERMS OF GRID POINT INDEX ON EACH 
Y-COORD INAIE LINE AT UP TO 5 STREAMWISE LOCATIONS 

IC 

CCOM 

INDEX REFERRING TO SUPERDIAGONAL HATRIX ELEMENTS 

ICDC(NN+1 ,2) 

CDC 

RECORD INDEX FOR READMS AND WRIIEMS MASS 
STORAGE DEV ICES-CDC COMPUTER ONLY 

ICOMP 

REST 

FLAG FOR COMPUTER OPTIONS 

ICONS ( 3 f NEQS ) 

OPER 

FLAG FOR CONVECTIVE FORMULATION BASED ON MACH 
NUMBER AND EQUATION 
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FORTRAN 

COHMON 

DESCRIPTION 

SYHBOL 

BLOCK 


ICORD 

GEOM 

FLAG FOR COORDINATE TRANSFORMATION OPTIONS 

ICPLB(NEQS f 2) 

EQN 

COUPLED EQUATION SENTINEL 

ID 

CCOM 

INDEX REFERRING TO VECTOR MATRIX ELEMENTS 

IDIFMO) 

OPER 

INDEX FOR TYPE OF DIFFERENCING OF BOUNDARY 
CONDITIONS 

IDMPY 

DMP 

DUMP LINE NO. IN Y DIRECTION 

IDMPZ 

DMP 

DUMP LINE NO. IN Z DIRECTION 

IDUM2(NN> 

ADDRG 

INDICATOR OF THE POINT IN DIFFERENCE MOLECULE 
WHERE Y DERIVATIVE IS TAKEN 

IDUH3(NN) 

ADDRG 

INDICATOR OF THE POINT IN DIFFERENCE MOLECULE 
WHERE Z DERIVATIVE IS TAKEN 

IEQ 

EQN 

EQUATION NUMBER INDEX 

IEQBC(20, 2,NEQS) 

BOUND 

SPECIFIED BOUNDARY CONDITION FOR WALL 

IEQNUM(NEQS) 

EQN 

IDENTIFICATION NUMBER OF EQUATIONS TO BE SOLVED 

IFBW(NN,NN> 

GRID 

INDICATOR WHICH TYPES GRID POINTS IN THE 
COMPUTATIONAL CROSS-SECTION PLANE. 

IFLARE 

OPER 

SENTINEL WHICH TELLS IF FLARE OPTION IS USED 

IGDMP 

DMP 

FLAG FOR DUMP OPTIONS 

IGEOM 

GEOM 

FLAG FOR COORDINATE SYSTEM OPTIONS 

IH 

GRID 

UPPER BOUNDARY POINT ON THE LINE WHERE IMPLICIT 
SOLUTION IS OBTAINED 

IHSTAG 

EQN 

FLAG FOR ENERGY EQUATION OPTIONS 

IJ 

DIFCOM 

COLUMN OR ROW NUMBER ON WHICH CALCULATION IS BEING 
MADE 

IL 

GRID 

LOWER BOUNDARY POINT ON THE LINE WHERE IMPLICIT 
SOLUTION IS OBTAINED 
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FORTRAN 

COHMON 

SYMBOL 

BLOCK 

IHIXL 

LAWU 

IND 

VAR 

INBC(NN) 

DIFCOM 

INDL 

VAR 

INH12 

ADD 

INH21 

ADD 

INH31 

ADD 

INH32 

ADD 

IOPTWF 

LAWW 

I0PIYZ(3,NEQS,2) 

OPER 

IPLOI 

CIO 

IPRINT 

CIO 

IPROF 

PRFILE 

IPRTE 

GASL 

IRSTIN 

REST 

IRSTOT 

REST 

ISONIC 

SUPER 


DESCRIPTION 

FLAG FOR MIXING LENGTH OPTIONS 

FLAG WHICH TELLS IF MESH POINT BELONGS TO 

SUPERSONIC OR SUBSONIC REGION 

MACH NUMBER INDICATOR 

STORAGE OF MACH NUMBER INDICATOR IN 

GENERAL PURPOSE STORAGE 

SENTINEL USED TO DETERMINE FORMULATION USED 
IN CALCULATION OF TRANSVERSE DERIVATIVE OF HI 
SENTINEL USED 10 DETERMINE FORMULATION USED 
IN CALCULATION OF STREAMWISE DERIVATIVE OF H2 
SENTINEL USED TO DETERMINE FORMULAI ION USED 
IN CALCULATION OF STREAMWISE DERIVATIVE OF H3 
SENTINEL USED TO DETERMINE FORMULATION USED 
IN CALCULATION OF TRANSVERSE DERIVATIVE OF H3 
SENTINEL WHICH DETERMINES WALL FUNCTION FORMULATION 
FLAG FOR DIFFERENCING FORMULATION BASED ON 
MACH NUMBER, EQUATION, AND ADI DIRECTION 
MARCHING STATION INTERVAL FOR STORAGE OF 
PLOTTING INFORMATION 

PRIMARY MARCHING STATION INTERVAL FOR PRINTING 

FLAG FOR INITIAL PROEILE OPTIONS 

FLAG WHICH DETERMINES EQ. OF STATE EORMULATION 

STREAMWISE STATION NUMBER FOR RESTART 

STREAMWISE INTERVAL FOR SAVING RESTART 

INFORMATION 

FLAG FOR SONIC LINE INTERPOLATION LOGIC 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


ISS(NN,5,4> 

SUPER 

GRID POINT LOCATION OF SONIC LINE 

ISSHFT 

SUPER 

INDEX USED 10 DETERMINE IF SONIC LINE IS LAST 
SUBSONIC POINI OR FIRSI SUPERSONIC POINT 

ISTART 

METRIC 

INITIAL CONDITION INDEX 

ISURF (4) 

BOUND 

LOCATION OF FOUR BOUNDARIES AS IDENTIFIED BY 
NSURF IN TERHS OF GRID POINT INDEX. ISURF(l) 

AND ISURF(2) ARE Y-INDICES OF NSURF(l) AND 
NSURF(2) RESPECT IVELY , WHEREAS ISURF<3) AND 
ISURFI4) ARE Z-INDICES OF NSURF<3) AND NSURF(4) 
RESPECTIVELY. 

ISW 

UNIVAC 

SENTINEL FOR WORD ADDRESSABLE OR SECTOR- 
ORIENTED MASS STORAGE DEVICE - UNIVAC ONLY 

ISYH 

ZPLOT 

SYMMEIRY OPTION FOR PLOTS 

ITRANS 

LAUW 

FLAG WHICH TELLS WHETHER TRANSITION TURBULENCE 
MODEL LOGIC IS USED 

IUNIIS 

UNITS 

FLAG USED TO DETERMINE SET OF DIMENSION UNITS USED 

IVARNO(NEQS) 

EQN 

IDENTIFICATION NUMBER OF DEPENDENT VARIABLES 

IVARPR ( 25) 

PRNT 

INDEX OF VARIABLES TO BE PRINTED 

IVISC 

vise 

FLAG FOR VISCOSIIY OPTIONS 

IUALF 

vise 

SENIINEl WHICH DETERMINES IF WALL FUNCTION 
LOGIC IS NEEDED IN THE CALCULATION OF WALL 
VISCOSITY 

IWR 

CIO 

SENIINEL FOR NAMELIST REST PRINT 

IYGD(NDIFM) 

ADDR 

GEOMETRY ADDRESSES 

IZCI<3) 

ZEX1 

RELATIVE UNIT NO. FOR VIRTUAL MEMORY STORAGE 

11 IG(8 ,3 ) 

FOCOM 

INDEX NEEDED IN THE CALCULATION OF METRIC 
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FORTRAN COMMON DESCRIPTION 

SYMBOL BLOCK 

INFORMATION 

12 IG(8 ,3) FGCOM INDEX NEEDE IN THE CALCULATION OF METRIC 

INFORMATION 

JA(5) ADDRF SHIFI LOGIC INDEX 

JADDO(NDIFM) ADDR ADDRESSES IN THE OPPOSITE DIRECTION 

JADDP(NDIFM) ADDR ADDRESSES IN THE PRIMARY DIRECTION 

JB0UND(NN,5,4) BOUND BOUNDARY IYPE INDICATOR AI EACH POINT ON BOUNDARY 

JDMAX EQN Y GRID POINT LOCATION OF MAXIHIUN DELTA 

JDRUM CIO ADD CODE DEVICE 

JDUM DIFCOM INDEX DENOTING RELATIVE POINT ABOUT WHICH 

DERIVATIVE IS LOCATED 

JEQBC(20,2,NEQS) BOUND SPECIFIED BOUNDARY CONDITION FOR NON-UALL 

JEQN(NEQS,2) EQN EXTERNAL EQUATION NUMBER 

JGSTOR LIN VALUE OF JG NEEDED BY SUBROUTINE EOS 

JMIN(LP,NN) BOUND LOWER GRID POINT INDEX LIMIT OF EACH COMPUTATIONAL 

LOOP(LP) ON EACH OF NN Y-COORDINATE LINES 
JPLOT CIO DEVICE FOR PLOTTING 

JPRINT CIO SECONDARY MARCHING STATION INIERVAL FOR PRINTING 

JPR0F(4) PRFILE SENT INEL FOR BOUNDARY VALUES DURING INITIAL PROFILE 

GENERATION 

JRSTIN REST LOGICAL FILE FROM WHICH RESIART INFORMATION 

IS READ 

JRSTOT REST LOGICAL F ILE ON WHICH RESTART INFORMATION 

IS WRITTEN 

JVAR(NEQS f 2) EQN VARIABLE NUMBER ASSOCIATED WITH EACH EQUATION 

DURING AN ADI SWEEP 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


JUR (5) 

CIO 

SENTINEL FOR NAMELIST PRINT OPTION 

JX 

OPER 

RELATIVE MARCHING STATION COUNTER 

JXDUM 

OPER 

ABSOLUTE MARCHING STATION COUNTER 

JXDUMP 

BHP 

MARCHING STATION WHEN DUMP OUTPUT IS REQUESTED 

KA(5) 

ABBRG 

INDICIES NECESSARY TO CALCULATE GEOMETRIC GROUPINGS 

KDMAX 

EQN 

Z GRID POINT LOCATION OF MAXIMIUM DELTA 

KBRUM 

CIO 

DEVICE FOR TEMPORARY STORAGE - USED IN GEOMETRY 
GENERATION 

KHAX(LP,NN) 

BOUND 

UPPER GRID POINT INDEX LIMIT OF EACH COMPUTATIONAL 
LOOP(LP) ON EACH OF NN Z-COORDINATE LINES 

KMIN(LP,NN> 

BOUND 

LOWER GRID POINT INDEX LIMIT OF EACH COMPUTATIONAL 
LOOP(LP) ON EACH OF NN Z-COORDINATE LINES 

JMAX(LP.NN) 

BOUND 

UPPER GRID POINT INDEX LIMIT OF EACH COMPUTATIONAL 
LOOP(LP) ON EACH OF NN Y-COORBINATE LINES 

LADD(3) 

BOUND 

ADDRESSES FOR BOUNDARY CONDITIONS 

LBRKY(IO) 

BOUND 

BREAK-OFF POINTS ON THE SPANWISE(Z-DIRECTION) 
COMPUTATIONAL BOUNDARY TO BE USED TO DETERMINE 
THE SURFACE NUMBER OF Y-PERSPECTIVE SEGMENTED 
BOUNDARY IN SUBROUTINE WHERE. 

LBRKZ(IO) 

BOUND 

BREAK-OFF POINIS ON THE IRANSVERSE( Y-DIRECTION) 
COMPUTATIONAL BOUNDARY TO BE USED TO DETERMINE 
THE SURFACE NUMBER OF Z-PERSPECI IVE SEGMENTED 
BOUNDARY IN SUBROUTINE WHERE. 

LBRUM 

CIO 

DEVICE FOR FINAL METRIC INFORMATION 

LEQ1 

EQN 

LOWEST INDEX OF EQUATIONS SOLVED EITHER BY 
COUPLED OR UNCOUPLED ADI SWEEP 
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FQRIRAN COMMON 

SYMBOL BLOCK 

LEQ2 EQN 

LEV(3) ADDR 

LEVEL ADDRG 

LEXIY(2 f NN) BOUND 


LEXTZ(2,NN) BOUND 


LGAKNN) ADDRG 

LGA2(NDIFM,NN) ADDRG 

LGA3(NBIFM,NN) ADDRG 

LGA4(NDIFM**2,NN) ADDRG 

LGA5(NDIFMAA2,NN) * ADDRG 


LP 

SECT 

LREF 

REF 

LREFI 

REF 

LSHFT 

GEOM 


DESCRIPTION 

HIGHEST INDEX OF EQUATIONS SOLVED EITHER BY 
COUPLED OR UNCOUPLED ADI SWEEP 
GEOMETRY LEVEL 
GEOMETRY LEVEL 

EXTENTS OF COMPUTATIONAL DOMAIN OF EACH 

Y-COORDINATE LINE 

LEXTY ( 1 ,NN) = LOWER LIMIT; 

LEXIY(2 f NN) = UPPER LIMII 

EXTENTS OF COMPUTATIONAL DOMAIN OF EACH 

Z-COORDINATE LINE 

LEXIZ(1,NN> = LOWER LIMIT; 

LEXTZ(2,NN) = UPPER LIMIT 

ADDRESS FOR POINT LOGIC OF GEOMETRIC VARIABLES 
ADDRESS FOR Y - DERIVATIVE OF GEOMETRIC 
VARIABLES 

ADDRESS FOR Z - DERIVATIVE LOGIC OF GEOMETRIC 
VARIABLES 

ADDRESS FOR CROSS DERIVATIVE(Y-Z) LOGIC OF 
GEOMETRIC VARIABLES 

ADDRESS FOR CROSS DERIVAIIVE(Z-Y) LOGIC OF 
GEOMETRIC VARIABLES 

INDEX OF COMPUTATIONAL LOOP IN EITHER Y OR Z 
DIRECTION. 

REFERENCE LENGTH 

INVERSE OE REFERENCE LENGTH 

SHIFT INDEX FOR COORDINATE TRANSFORMATION 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


LVG(8) 

BOUND 

BOUNDARY POINI INDICATOR FOR BOUNDARY CONDITION 

MASS1 

CIO 

GENERAL PURPOSE MASS STORAGE DEVICE 

MASS2 

CIO 

GENERAL PURPOSE MASS STORAGE DEVICE 

MCPLD 

EQN 

NUMBER OF COUPLED EQUATIONS TO BE SOLVED 

HEFF 

VAR 

INDEX FOR EFFECTIVE VISCOSIIY 

MEFFL 

VAR 

STORAGE LEVEL OF EFFECTIVE VISCOSITY IN 
GENERAL PURPOSE STORAGE 

MEQK 

EQN 

LEQ1 - 1 

MEQS 

EQN 

TOTAL NUMBER OF EQUATIONS TO BE SOLVED 

MEQS1 

EQN 

INDEX OF FIRST EQUATION TO BE SOLVED 

MEQS2 

EQN 

INDEX OF LAST EQUATION TO BE SOLVED 

MQDMP 

DMP 

FLAG FOR DUMP OPTIONS 

MGD1 

GRID 

IL + 1 

MGD2 

GRID 

IH - 1 

MIN 

CIO 

INPUI DEVICE 

MINE 

FREE 

FREE STREAM HACH NUMBER 

ML 

VAR 

INDEX FOR MIXING LENGTH 

HLEVEL 

PARAM 

MAXIMUM NO. OF STORAGE LEVELS 

MLL 

VAR 

STORAGE LEVEL OF MIXING LENGTH IN 
GENERAL PURPOSE STORAGE 

MN 

VAR 

INDEX FOR MACH NUMBER 

MNL 

VAR 

STORAGE LEVEL OF MACH NO IN 
GENERAL PURPOSE STORAGE 

MOUT 

CIO 

OUTPUT DEVICE 

MREF 

REF 

REFERENCE MACH NUMBER 

HREF I 

REF 

INVERSE OF REFERENCE MACH NUMBER 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

MSDD 

CIO 

TEMPORARY MASS STORAGE DEVICE 

MSD1 

CIOD 

GENERAL PURPOSE MASS STORAGE DEVICE 

MSD2 

CIOD 

GENERAL PURPOSE MASS STORAGE DEVICE 

MSECY 

SECT 

MAXIMUM NO. OF Y-PERSPECTIVE SUBSECTIONS 
(USED 10 SEI UP THE SUBSECTIONS FOR CROSEC) 

MSECZ 

SECT 

MAXIMUM NUMBER OF Z-PERSPECTIVE SUBSECTIONS 
(USED 10 SET UP THE SUBSECTIONS FOR CROSEC) 

MSGVAR(25) 

PRNI 

TITLE OF VARIABLES TO BE PRINTED 

MU 

VAR 

INDEX FOR LAMINAR VISCOSITY 

MUINF 

FREE 

FREE STREAM LAMINAR VISCOSITY 

MUL 

VAR 

STORAGE LEVEL OF LAMINAR VISCOSITY IN 
GENERAL PURPOSE STORAGE 

MUREF 

REF 

REFERENCE VISCOSIIY 

MUREFI 

REF 

INVERSE OF REFERENCE VISCOSITY 

MUT 

VAR 

INDEX FOR TURBULENT VISCOSITY 

MUIL 

VAR 

STORAGE LEVEL OF TURBULENT VISCOSIIY IN 
GENERAL PURPOSE STORAGE 

MW INF 

FREE 

FREE STREAM MOLECULAR WEIGHT 

MWREF 

REF 

REFERENCE MOLECULAR WEIGHT 

MUREFI 

REF 

INVERSE OF REFERENCE MOLECULAR WEIGHT 

MZVAR 

PARAM 

MAXIMUM NUMBER OF STORAGE VARIABLES 

NABC 

CCOM 

NABC = ID 

NANG 

METRIC 

ANGLE OF COORDINATE LINES RELATIVE TO HORIZONTAL 

NBRKX 

BOUND 

NUMBER OF STREAMU ISE LOCATIONS WHERE SOLID BODIES 
EMBEDDED IN THE COMPUTATIONAL DOMAIN START OR 
TERMINATE. 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


NCMAX( 5) 

SECT 

UPPER LIMIT OF BOUNDARY LINE WHERE ENDCAP 
SOLUTION IS OBTAINED. 

NCHIN(S) 

SECT 

LOWER LIMIT OF BOUNDARY LINE WHERE ENDCAP 
SOLUTION IS OBTAINED. 

NCPLD 

PARAH 

HAXIMIUN NUMBER OF COUPLED EQUATIONS 

NCPLD2 

CCOM 

NCPLDAA2 

NCRF 

SECI 

EQUAL TO EITHER NYCRF(ISEC) OR NZCRF(ISEC) WHERE 
ISEC INDICATES A SUBSECTION IN THE CROSS-SECTION 

NCRL 

SECT 

EQUAL TO EITHER NYCRL(ISEC) OR NZCRL(ISEC) 
WHERE ISEC INDICATES A SUBSECTION IN THE 
CROSS-SECTION. 

NCTR 

BIFCOM 

CENTER OE DIFFERENCE MOLECULE 

NCUP 

EQN 

NUMBER OF COUPLED EQUATIONS 

NDIFM 

PARAM 

MAXIMIUM NUMBER OF GRID POINTS IN A DIFFERENCE 
MOLECULE 

NBIFMT 

BIFCOM 

2ANDIFM 

NDIFM1 

BIFCOM 

NDIFM - 1 

NDIFP1 

BIFCOM 

NDIFH + 1 

NE(2) 

GEOH 

NUMBER OF GRID POINTS IN Y AND Z DIRECTIONS 

NEQN<NEQS,2> 

EQN 

NUMBER OF COUPLED EQUATIONS TO BE SOLVED IN EACH 
ADI SWEEP 

NEQS 

PARAM 

MAXIMIUM NUMBER OF EQUATIONS IN CODE 

NEY 

GRIB 

NUMBER OF GRID POINTS IN Y DIRECTION 

NEYM1 

GRIB 

NEY - 1 

NEZ 

GRIB 

NUMBER OF GRID POINTS IN Z DIRECTION 

NEZM1 

GRID 

NEZ - 1 
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FORTRAN 

SYHBOL 

NE2S 

NFILE 

NGEOMV 

NH1 

NH12 

NH2 

NH21 

NH3 

NH31 

NH32 

N I IT 

NIN 

NINC 

NJD 

NJF(IO) 

NJL(IO) 

NJSEG 

NKF(IO) 

NKL(IO) 

NKSEG 


COMMON DESCRIPTION 

BLOCK 

CIO SENTINEL FOR SPREADING OF 2-D PROFILE TO 3-D 

REST SEQUENCE NUMBER OF RESTART INFORMATION 

METRIC NUMBER OF HEIRIC COEFFICIENTS AND DERIVATIVES 
MEIRIC INDEX FOR METRIC COEFFICIENT IN X-DIRECIION 
METRIC INDEX FOR DERIVATIVE OF X METRIC IN Y DIRECTION 
METRIC INDEX FOR METRIC COEFFICIENT IN Y-DIRECTION 
METRIC INDEX FOR DERIVATIVE OF Y METRIC IN X DIRECTION 
MEIRIC INDEX FOR METRIC COEFFICIENT IN Z-DIRECTION 
MEIRIC INDEX FOR DERIVATIVE OF Z METRIC IN X DIRECTION 

METRIC INDEX FOR DERIVATIVE OF Z METRIC IN Y DIRECTION 

CIO NO. OF FALSE MARCHING STEPS USED TO GENERATE 

THE INITIAL PROFILE 

DIFCOM MAXIMUM LINES OF STORAGE IN CORE AT ONE TIME. 
DIFCOM NIN/2 + I 

DIFCOM GRID POINT LOCATION FOR START OF SECOND SWEEP 

SECT LEFT BOUNDARY INDEX OF EACH Y-PERSPECT IVE 

SUBSECTION USED FOR Y-SWEEP OF ADI. 

SECI RIGHT BOUNDARY INDEX OF EACH Y-PERSPECT IVE 

SUBSECTION USED FOR Y-SWEEP OF ADI. 

SECT NUMBER OF SUBSECTIONS IN EACH CROSS SECTION 

USED FOR Y-SWEEP OF ADI 

SECT LOWER BOUNDARY INDEX OF EACH Z-PERSPECT IVE 

SUBSECT ION USED FOR Z-SWEEP OF ADI. 

SECT UPPER BOUNDARY INDEX OF EACH Z-PERSPECT IVE 

SUBSECTION USED FOR Z-SWEEP OF ADI. 

SECT NUMBER OF SUBSECTIONS IN EACH CROSS SECTION 


78 


FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

USED FOR Z-SUEEP OF ADI. 

NLOOPY(NN) 

SECT 

NUMBER OF LOOPS ON Y-COORDINAIE LINE 

NLOOPZ(NN) 

SECT 

NUMBER OF LOOPS ON Z-COORBINATE LINE. 

NHAXUF 

LAUU 

MAXIMIUM NUMBER OE ITERATIONS ALLOWABLE IN 
CALCULAIIOON OE WALL SHEAR VELOCIIY 

NN 

PARAM 

MAXIMIUM NUMBER OF GRID POINTS IN Y OR Z DIRECTION 

NPADI 

EQN 

NUMBER OF COUPLED AND UNCOUPLED EQUATIONS TO BE 
SOLVED DURING EACH ADI SWEEP 

NPROFUO) 

PRFILE 

NUMBER OF SURFACES WHERE INITIAL BOUNDARY LAYER 
PROFILES ARE GENERATED IN EACH LOOP EITHER IN 
Y OR Z DIRECTION. NPROF IS EITHER 1 OR 2. 

NPTSX(4) 

INTBC 

NUMBER OF SIREAMUISE LOCATIONS WHERE TRANSPIRATION 
DATA IS INPUI 

NPTSYZ( 15,4) 

INTBC 

NUMBER OF CROSS-PLANE LOCATIONS UHERE TRANSPIRATION 
IS INPUI 

NPUNCH 

CIO 

PUNCH DEVICE 

NRGT(2) 

DIFCOM 

NCTR POINTS FROM RIGHT OR TOP BOUNDARY 

NS 

GEOM 

LAST MARCHING STATION 

NSAVED 

RESI 

SEQUENCE NUMBER OF RESTART STATIONS SAVED 

NSECRY 

SECT 

NUMBER OF Y-PERSPECT IVE SUBSECTIONS IN THE 
CROSS-SECTION FOR SUBROUTINE CROSEC. 

NSECRZ 

SECT 

NUMBER OF Z-PERSPECT IVE SUBSECTIONS IN THE 
CROSS-SECTION FOR SUBROUTINE CROSEC. 

NSMAX 

PARAM 

2 GREATER THAN NS 

NSURF (4) 

BOUND 

SURFACE NUMBER OF FOUR BOUNDARIES WITH MINIMUM 
DISTANCE AWAY FROM A GRID POINT CONSIDERED. 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 

NSURF(l) AND NSURF(2) INDICATE Y-PERSPECT IVE 
SURFACE NUMBER r UHEREAS NSURF(3) AND NSURF<4) 
INDICATE Z-PERSPECTIVE SURFACE NUMBER 

NUNERR 

CIO 

DEVICE FOR MGAUSS ERROR CHECK 

NVSOLV 

EQN 

NUMBER OF DEPENDENT VARIABLES 

NWORD2(50) 

STRAGE 

SIZE OF EACH COMMON BLOCK 

NYCRF(IO) 

SECT 

LOWER BOUNDARY INDEX OF EACH Z-PERSPECTIVE 
SUBSECTION FOR USE IN CROSEC. 

NYCRLUO) 

SECT 

UPPER BOUNDARY INDEX OF EACH Z-PERSPECTIVE 
SUBSECTION FOR USE IN CROSEC. 

NZCRF(IO) 

SECI 

LEFT BOUNDARY INDEX OF EACH Y-PERSPECTIVE 
SUBSECTION FOR USE IN CROSEC. 

NZCRL(IO) 

SECT 

RIGHT BOUNDARY INDEX OF EACH Y-PERSPECIIVE 
SUBSECTION FOR USE IN CROSEC. 

OMBWB 

LIN 

1.0 - BUD 

OMEGWF 

LAWW 

UNDER-RELAXATION FACTOR FOR WALL FUNCTION 
FORMULATION 

P 

VAR 

INDEX FOR STATIC PRESSURE 

PC0N1 

REF 

(GAMMA-1. 0)/GAMMA 

PC0N2 

REF 

0.5 A PC0N1 

PINF 

FREE 

FREE STREAM STATIC PRESSURE 

PL 

VAR 

STORAGE LEVEL OF STATIC PRESSURE IN 
GENERAL PURPOSE STORAGE 

PLTFLD(NN,NN f 8) 

ZPLOT 

GENERAL PURPOSE STORAGE FOR PLOT INFORMATION 

PR 

REF 

PRANDTL NUMBER 

PREF 

REF 

REFERENCE PRESSURE 
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FORTRAN 

SYMBOL 

COMMON 

BLOCK 

DESCRIPTION 

PREFI 

REF 

INVERSE OF REFERENCE PRESSURE 

PREPS 

vise 

PRANDTL NO. IN TURBULENT ENERGY DISSIPATION 



EQUATION 

PRESS(20,2) 

BOUND 

SPECIFIED PRESSURE AT BOUNDARY 

PRI 

REF 

TURBULENT PRANDTL NO. 

PRIKE 

vise 

PRANDTL NO. IN TURBOLENI KINECIIC ENERGY EQUATION 

PZERO 

FREE 

STAGNATION PRESSURE 

Q1D(8,NN> 

METRIC 

INTERMEDIATE STORAGE ARRAY FOR METRIC 



INFORMATION 

Q2D(8,NN) 

METRIC 

INTERMEDIATE STORAGE ARRAY FOR METRIC 



INFORMATION 

R 

VAR 

INDEX FOR DENSITY 

RATLD 

LAWU 

EMPIRICAL NUMERICAL CONSTANT IN MIXING LENGTH 



COMPUTATION 

RE 

REF 

REYNOLDS NUMBER 

REI 

REF 

INVERSE OF REYNOLDS NUMBER 

REI2 

REF 

2.0 * REI 

REPL 

FREE 

REYNOLDS NUMBER PER UNIT LENGTH 

RGAS 

REF 

GAS CONSTANT 

RHO(NDIFM) 

ADBR 

STORAGE OF DENSITY FOR EIRST DERIVATIVES 

RHOINF 

FREE 

FREE STREAM DENSITY 

RHOREE 

REF 

REFERENCE DENSITY 

RH0UL(NN,5,4> 

BOUND 

STORAGE FOR DENSITY ON THE BOUNDARY 

RHREFI 

REF 

INVERSE OF REFERENCE DENSITY 

RUN IV 

UNITS 

UNIVERSAL GAS CONSTANT 

SAVE(NEQS,NN> 

EQN 

STORAGE FOR CHANGES DURING FIRST ADI SWEEP 
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FORIRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


SN(NN) 

LIN 

STORAGE FOR N IH LEVEL TERMS 

SQ1 

ADD 

ADD CODE STREAMWISE LOCATION 

SQ2 

ADD 

ADD CODE STREAMUISE LOCATION 

SYSTEM 

ZPLOT 

SENTINEL FOR COORDINATE SYSTEM - PLOTS ONLY 

T 

VAR 

INDEX FOR STATIC TEMPERATURE 

TEMPS(NN f 5,4) 

BOUND 

STORAGE ARRAY FOR TEMPERATURE ON BOUNDARY 

IEMPSN(NN f 5 f 4) 

BOUND 

STORAGE ARRAY FOR TEMPERATURE AT NTH STREAMUISE 
STATION - BOUNDARIES ONLY 

TINE 

FREE 

FREE STREAM STATIC TEMPERATURE 

TITLE<6) 

ZPLOT 

TIILE FOR PLOT FILE 

TKE 

VAR 

INDEX FOR TURBULENI KINEIIC ENERGY 

TKEINF 

LAUU 

FREE STREAM TURBULENT KINETIC ENERGY 

IL 

VAR 

STORAGE LEVEL OF STATIC TEMPERATURE IN 
GENERAL PURPOSE STORAGE 

IREF 

REF 

REFERENCE IEMPERATURE 

TREE I 

REF 

INVERSE OF REFERENCE TEMPERATURE 

TT1 (2) 

GEOM 

MESH DISTRIBUTION FACTOR 

TT2(2) 

GEOH 

MESH DISTRIBUTION EACTOR 

TUALL<20,2) 

BOUND 

SPECIFIED TEMPERATURE AI BOUNDARY 

TUOD 

GEOH 

FLAG FOR TWO DIMENSIONAL LOGIC 

TZERO 

FREE 

STAGNATION TEMPERATURE 

T2(2) 

GEOM 

MESH DISTRIBUTION EACTOR 

U 

VAR 

INDEX FOR VELOCITY IN X-DIRECTION 

UDUE(20 r 2) 

LAUU 

BOUNDARY LAYER THICKNESS SAMPLING CRIIERIA 

UDUM(NBIFM) 

ADDR 

STORAGE OF VELOCITY FOR FIRST DERIVATIVES 

UINF 

FREE 

FREE STREAM VELOCITY 
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FORIRAN 

SYMBOL 

UREF 

UREFI 

USCALE 

USTAR(NN,5,4> 
UT IL(NDIFM) 

V 

VELSQ(NN,NN) 

VKB 

VKC 

VNO( 15,15,4) 

U 

X(502) 

XBRKX(S) 

XENIR 

XG1(NN,2) 

XG2(NN, 2) 

XOB ( 4 ) 

XVNO (15,15,4) 
XO(2 ) 

Y(NN,NN,2) 


COHMON DESCRIPTION 

BLOCK 

REF REFERENCE VELOCIIY 

REF INVERSE OF REFERENCE VELOCITY 

ADD METRIC SCALE FACIOR 

LAUU FRICTION VELOCITY ON SOLID UALL BOUNDARY 

ADDR STORAGE OF DEPENDENT VARIABLE FOR FIRST DERIVATIVES 

VAR INDEX FOR VELOCIIY IN Y-DIRECIION 

STORAGE FOR UAA2 + VAA2 + UAA2 - EQUIVALENCED TO 
PLTFLD< 1,1,4) 

LAUU SECOND CONSTANT IN LOGARITHMIC LAU OF THE UALL 

LAUU VON KARMAN CONSTANT 

INTBC INPUT TRANSPIRATION RATES 

VAR INDEX FOR VELOCIIY IN Z-DIRECTION 

GEOM STREAMUISE LOCATION 

BOUND STREAMUISE LOCATION UHERE AN EMBEDDED SOLID 

BODY STARTS OR TERMINATES 
GEOM STARTING STREAMUISE LOCATION 

OPER FIRST DERIVATIVES OF COMPUTATIONAL COORDINATES 

UITH RESPECT TO PHYSICAL COORDINATES 
OPER SECOND DERIVATIVES OF COMPUTATIONAL 

COORDINATES UITH RESPECI TO PHYSICAL COORDINATES 
BOUND INITIAL LOCATION OF SOLID OBSTACLE IN X DIRECTION 
INTBC STREAMUISE LOCATIONS UHERE TRANSPIRATION DATA IS 
INPUT 

GEOM MESH DISTRIBUTION FACTOR 

PHYSICAL DISTANCES FROM BOUNDARIES - EQUIVALENCED 
TO PLTFLD < 1 f 1 , 2 ) 
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FORTRAN 

COMMON 

DESCRIPTION 

SYMBOL 

BLOCK 


YAW 

REE 

YAW ANGLE 

YPLUSL (NN f NN) 


STORAGE FOR RHO * UIAU / VISLAM - EQUIVALENCED TO 
PLTFLD( 1 ,1,1) 

YS(2,2) 

GEOH 

NONDIMENS IONAL EXTENTS OF COMPUTATIONAL DOMAIN 

YSAVE(NN,2) 

GEOM 

COMPUTATIONAL COORDINATES 

YZPROF(NN) 

PRFILE 

TEMPORARY STORAGE ARRAY FOR PHYSICAL COORDINATES 

YZUN0(15,15,4) 

INTBC 

CROSS-PLANE LOCATIONS WHERE IRANSPIRAIION DATA IS 
INPUT 

ZNIRN<NDIM,3> 

ZEX1 

ARRAY FOR VIRTUAL MEMORY STORAGE 
NDIM = NN A NN A MZVAR A MLEVEL 

ZZ(M,L,K) 

VARZZ 

GENERAL PURPOSE STORAGE FOR DEPENDENT AND DERIVED 
VARIABLES - M = MZVAR, L = MLEVEL, K = NDIFM A NN 
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Sample Input and Sample Output 


Tables I and II present sample input and output for one of the test 
cases previously discussed. The sample input was used to calculate the flow 
field of rectangular supersonic jet (aspect ratio = 1) interacting with the 
ambient stream. The input data is for an initial run (IRSTIN = 0) with a 
restart to be written every 10 marching steps (IRSTOT=10). The case is to be 
run on CRAY-1 computer. The streamwise and two transverse momentum as well 
as continuity and stagnation enthalpy (IHSTAG =1) version of the energy 
equation are to be solved. The reference length is 2.0m (IUNITS * 2), the 
Reynold's number per m is 1.6685 x 10 6 , the free stream Mach number 2.0 and 
the free stream pressure is 2864.0 Nt/m . 50 x 50 grid points are utilized 

in both directions of the computational cross-section (NE = 50,50) and a 
Cartesian coordinate system is to be utilized ( IGECM = 1). The initial run 
is to be marched 10 steps (NS = 10) starting at a streamwise location of 0.0 
(XENTR « 0.0). The streamwise marching step size is 0.01 (DELX = 0.01). 

Then grid point characteristics is specified (IFBW(1,1) =....) to define 
a computational domain in the cross-section. For a detailed explanation of 
the indicators, computer program list should be consulted. During this 
initial run grid point characteristics will change at one streamwise location 
(NBRKX - 1). At the streamwise location 0.015 (XBRKX(l) = 0.015) an embedded 
solid body will begin (IED GE = 0). The solid body is a splitter bounded by 
two surface-fitted coordinate lines (18 and 19 in Y-direction and 18 and 19 
in Z-direction) . In the Namelist LISTR, IFBW is respecified as shown in 
Table I. An information on IBOUND and JEQBC corresponding to the new IFBW is 
also provided. Grid points are clustered in the vicinity of the embedded 
solid body, i.e., about 0.2 in Y-direction and 0.2 in Z-direction (X0(1) * 
0.2, 0.2). The packing is to be moderately tight (T2(l) - 3.0, 3.0). An 
initial profile is generated within the program to supply a uniform stream 
needed for this case (IPROF = 1) and modified to give two different 
magnitudes of velocity. Printout is given every 5 steps (IPRINT = 5) and 
plot information is written every 5 steps (IPL0T = 5). 

The output for the three-dimensional case consists of NAMELIST 
information and difference weight information and flow fieLd information at 
each 5th streamwise station. The NAMELIST information is provided as a means 
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for the user to check the input data. Three-dimensional flow field output is 
controlled by the variable IVARPR. For three dimensions, however, the output 
is in the form of a cross-sectional plane of output. All variables (except 
the pressure) are in a non-dimensional form with respect to the reference 
conditions which are displayed in NAMELIST LIST 2. The pressure terms are 
nondimensionalized with respect to the freestream pressure. The integer 
variables IY and IZ represent the transverse and spanwise grid point 
locations respectively while Z and Y are the corresponding computational 
positions. Table II is a portion of the output for the three-dimensional 
rectangualr jet case i.e., the cross plane distribution of the streamwise 
velocity (UVEL) and pressure (PRES). Other variables can (and were) printed 
out, but for reasons of economy of space are not presented here. Following 
the flow field information is the subsonic-supersonic grid point position 
indicator (ISS) with respect to the lower and upper surface (right and left 
surface for IADI - 3), and the boundary indicator ( JBOUND) . The ISS values 
tell the grid point where the flow transitions from subsonic to supersonic 
flow while the variable JBOUND tells the type of surface (JBOUND = 1 
corresponding to a wall and JBOUND = 2 corresponding to a nonwall). Finally, 
plot file information is displayed. 
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SUBSONIC EXTERNAL STEAM 



SUPERSONIC EXTERNAL STEAM 

PLUME BOW 



FIGURE 2 - SCHEMATIC OF JET NEAR FIELD STRUCTURE FOR 
SUPERSONIC AND SUBSONIC EXTERNAL FLOWS 
(REF. 18) 
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FIGURE 3 - SCHEMATIC OF AXISYMMETRIC JET FLOW 
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FIGURE 4(b) - STREAMWISE STATIC PRESSURE DISTRIBUTION 




MIXING REGION 


MIXING REGION 



SIDE VIEW 



FRONT VIEW 


FIGURE 5 - SCHEMATIC OF 3- DIMENSIONAL JET FLOW 
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FIGURE 6 - SCHEMATIC OF COMPUTATIONAL DOMAIN (SIDE VIEW) 
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FIGURE 7 - CROSS-SECTION OF 
(ASPECT RATIO = I) 
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TABLE I - Sample Three-Dimensional Input 


TOP RECORD 

N0ZZLE(AR = 1 ) 

1 1.0 

SREST 

ICOMP = 3, 

IRSIIN = 0, NFILE = 0, NSAOED = 0, 
IRSTOT = 10, 

JRSTOT ■ 10, 


SEND 
4L 1ST 1 
IHSTAG = 1, 
IBOUND = 40*2, 
IEQBC ( 1,1,1) = 
IEQBC (1,1,2) = 
IEQBC (1,1,3) = 
IEQBC (1,1,4) = 
IEQBC (1,1,5) = 
JEQBC (1,1,1) = 
JEQBC ( 1,1,2) = 
JEQBC (1,1, 3) = 
JEQBC ( 1,1,4) = 
JEQBC ( 1,1,5) = 
SEND 


40*2, 

40*2, 

40*2. 

40*16, 

40*8, 

11,12,38*11, 

2,12,18*2,20*11, 

11,12,18*11,20*2, 

14,12,38*14, 

17,12,38*17, 


&LIST2 
LREF = 1.0, 

REPL = 1 . 6685E+06 , 
MINF = 2.0, 

PINF = 2864.0, 

PR = 0.71, 

IUN ITS = 2, 

SEND 

&LIST3 

TWOD = .FALSE., 

X0< 1 ) = 0.2, 0.2, 

T2 ( 1 ) = 3. 0,3.0, 

NE = 50,50, 

IGEOM = 1, 

DELX = 0.01, 

IAP = 1, 

DXM IN = 0.01, 

DXMAX = 0.01, 


AP = 1.0, 

NS = 10, 
XENTR = 0.0 
IFBU (1,1) = 


5,48*3,5,50*2, 

3,48*1,3,50*2, 

3,48*1,3,50*2, 

3,48*1,3,50*2, 



3,48*1, 3,50*2 , 
3,48*1 ,3,50*2, 
3 , 48A1 , 3 , 50*2, 
3 , 48A1 , 3,50*2 , 
3 f 48A1 , 3,50*2 , 
3,48*1,3,50*2, 
3 , 48A1 f 3 r 50*2 , 
3 , 48A1 , 3,50*2 , 
3 , 48A1 , 3,50*2 , 
3 , 48A1 , 3, 50A2 , 
3 , 48A1 ,3,50*2, 
3 f 48A1 ,3,50*2, 
3 , 48A1 ,3,50*2, 
3,48*1, 3,50*2, 
3 , 48A1 , 3 , 50*2 , 
3 , 48A1 ,3,50*2, 
3 , 48A1 ,3,50*2, 
3,48*1,3,50*2, 
3, 48A1 ,3,50*2, 
3,48*1,3,50*2, 
3, 48A1 ,3,50*2, 
3 f 48A1 ,3,50*2 , 
3 , 48A1 , 3 , 50A2 , 
3,48*1 ,3,50*2 , 
3, 48*1,3,50*2, 
3,48*1 ,3,50*2 , 
3, 48A1 ,3,50*2, 
3 , 48A1 ,3,50*2-, 
3 , 48A1 ,3,50*2, 
3, 48A1 ,3,50*2 , 
3 , 48A1 , 3 , 50*2 , 
3 , 48A1 ,3, 50A2 f 
3 r 48A1 ,3,50*2, 
3 , 48A1 ,3, 50*2 , 
3 , 48A1 , 3 , 50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
3,48*1,3,50*2, 
5,48*3,5,50*2, 


NBRKX = 1, 

XBRKX<1) = 0.015, 

IBRKMN = 17*18,2*1,481*0, 

IBRKMX ■ 19*19,481*0, 

IEDQE = 0, 

SEND 
&LIST4 
IVISC = 3, 

IPROE = 1, 

IMIXL = 1, 

SEND 
&L IST5 
IPRINT » 5, 

IPLOT » 5, 

SEND 

SLISTR 

IFBW (1,1) = 5,16*3,2*5,30*3,5,50*2, 
3,16*1,2*3,30*1 ,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3 ,16*1,2*3,30*1,3,50*2, 
3,16*1, 2*3,30*1,3,50*2, 
3, 16*1 ,2*3,30*1 ,3, 50-*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 
3,16*1,2*3,30*1,3,50*2, 


3, 

5 , 

5, 

IBOUND = 2,2*1 
JEQBC <1,1,1) * 
JEQBC (1,1,2) = 
JEQBC (1,1,3) = 
JEQBC (1,1,4) = 
JEQBC (1,1,5) « 


16*1,2*3,30*1,3,50*2, 

16*3,5,3,30*1,3,50*2, 

17*3,4,30*1,3,50*2, 

,2,3*1,3*2,10*2,2,2*1,2,3*1,3*2,10*2 

3*11,12,3*11,12,11,12,30*11, 

3*2,12,3*2,12,2,12,10*2,20*11, 

3*11,12,3*11,12,11,12,10*11,20*2, 

3*14,12,3*14,12,14,12,30*14, 

3*17,12,3*17,12,17,12,30*17, 


T 


SEND 

EOF 
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